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FOREWORD 


The  work  described  in  this  report  was  performed  by  Dr.  Russell  H.  Lyddane  of  Ebeling 
Associates,  Inc.,  Scotia,  New  York.  Given  are  the  mathematical  derivations  and  the  complete 
FORTRAN  Code. 

The  Lyddane  program  is  estimated  to  operate  at  least  two  orders  of  magnitude  faster  than 
standard  orbit  integration  techniques.  It  is  eminently  suitable  for  the  computation  of  satellite 
orbits  with  moderate  accuracy  over  long  periods  of  time. 
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THOMAS  A.  CLARE,  Head 
Strategic  Systems  Department 
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INTRODUCTION 


The  General  Perturbation  Satellite  (GPS)  computer  program  was  designed  to  conduct  the 
integration  of  the  equations  of  motion  of  an  earth  satellite  under  the  influence  of  (1)  an 
arbitrary  zonal-harmonic  gravitational  field  of  the  earth  and  (2)  the  gravitational  attractions 
of  the  sun  and  moon.  The  zonal  harmonics  (as  the  program  stands)  go  up  to  the  twentieth, 
but  additional  terms  can  be  added  by  changing  a  few  parameters  at  the  start  of  the  program. 
There  is  no  provision  for  tesseral  harmonics.  Sun  and  moon  coordinates  as  functions  of  time 
are  read  from  a  data  file  that  will  not  be  described  here,  since  its  construction  is  outside 
the  program. 

The  program  employs  Poincare-like  canonical  variables  and  transforms  them  by  general 
perturbation  theory  (von  Zeipel)  to  remove  first-order  short-period  terms  (higher-order  short- 
period  terms  are  neglected).  The  equations  in  the  mean  variables  thus  produced  are  integrated 
numerically,  with  a  large  step  size  (set  currently  at  1/4  day).  Since  long-period  terms  receive 
no  special  treatment,  there  is  nothing  special  about  the  critical  angle;  and  since  the  variables 
are  Poincare-like,  there  are  no  difficulties  at  low  eccentricities  or  inclinations. 

A  listing  of  the  GPS  program,  in  FORTRAN  IV,  is  shown  in  Appendix  A,  with  line 
numbers  and  with  concordances  for  the  various  subdivisions.  This  is  referenced  as  appropriate 
to  relate  derivations  to  the  final  program.  The  relationships  among  the  subroutines  and  common 
variables  are  shown  in  Appendix  B. 


COORDINATES,  ELEMENTS,  AND  CANONICAL  VARIABLES 


The  basic  Cartesian  coordinate  system  employed  is  the  nonrotating  system  with  center 
at  the  earth’s  center,  Z-axis  along  the  north  polar  axis,  epoch  1950.0.  This  is  the  system  in 
which  solar  and  lunar  coordinates  are  also  given.  The  coordinates  of  the  satellite  are  x,  y,  z 
and  the  velocity  components  v#,  vy,  vf. 
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The  (osculating)  elliptic  elements  are  the  usual  set: 

a  =  semi-major  axis 

e  =  eccentricity 

i  .=  inclination 

M  =  mean  anomaly 

u  =  argument  of  pericenter 

ft  -  aigument  of  node 

The  auxiliary  variables  follow: 
f  =  true  anomaly 
E  =  eccentric  anomaly 
oj  =  cj  +  ft  argument  of  longitude 
a  -  E  +  to  eccentric  longitude 
u  =  f  +  u>  true  longitude 

The  Delaunay  canonical  variables,  used  only  as  convenient  intermediate  variables,  are  also 
as  usual: 

L  =  Vpa  C  =  M 

G  =  L  V 1  -  ei  g  =  w 

H  =  G  cos  i  h  = 

where  p  is  the  gravitational  constant.  Usually  8,  g,  and  h  will  be  employed  instead  of  M.  co. 
and  ft. 
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The  canonical  variables  used  in  the  program  follow: 


X  =  8  +  gj 


i  =  V2(L-G)  cos  to  7j  =  -  v  2(L-G)  sin  gj 

a  =  V 2(G-H)  cos  h  r  =  -  V  2(G-H)  sin  h 

Of  course,  X  is  the  mean  longitude.  Proof  that  this  is  indeed  a  canonical  transformation 
is  easily  derived  (Reference  1). 

Several  more  intermediate  quantities  are  convenient  to  use  with  these: 


f  =  G/L  =  v  1  -  e2  =  l  -  (^2  +  v2  )/2L 


ec  =  e  cos  gj 


■VW 


e„  =  e  sin  gj  = 


(-  7?) 


p  =  H/G  =  cos  i  =  1  -  ( a 2  +  t2  )/  2  f  L 


Ic  =  sin  i  cos  h  =  a  V ( 1  +  p)/2Lf 


ls  =  -  sin  i  sin  h  =  r  V ( 1  +  p/2L  f 


Details  of  the  transformations  between  these  sets  of  coordinates  are  given  in  Appendix  C. 


H.  Goldstein,  “Classical  Mechanics" 
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MEAN  ELEMENTS:  THE  VON  ZEIPEL  TRANSORM ATION 


The  perturbing  potential  in  the  Hamiltonian  is  dominated  by  the  contribution  of  the  second 
zonal  harmonic,  whose  coefficient  is  the  small  perturbing  parameter  (~  0.001).  All  other  terms 
in  the  perturbing  potential  are  of  the  order  of  the  square  of  this  parameter.  We  expand  the 
Hamiltonian  in  powers  of  the  parameter: 

F  =  1-  +  F  +  F  F  =  /li2/2L2 

0  12  0  ^  ' 

and  introduce  a  determining  function  S  expanded  in  the  same  way: 

S  =  S()  +  S,  +  S,  where  SQ  =  L'X  +  j-'r)  +  o't 

and  the  primed  quantities  are  the  mean  elements. 

We  will  choose  S  so  that  the  first-order  short-period  terms  are  transformed  out.  Higher- 
order  short-period  terms  will  be  neglected,  except  for  taking  into  account  their  effect  on  the 
initial  value  of  L*.  as  described  below. 


Since  F(L.  £.  o.  X.  77.  r)  =  F  ( L  .  S  ,  0  .  - .  77.  r ). 


,  as\  „  /as  as  as  , 

F„  (  —  )  +  F  ( .  —  .  —  .  A .  77 

0  \a,  /  1  \ax  377  a r 


•  +  f2  - 

.  ,  *  /  ,  ,  ,  as  as  \  * 

,  +  I  .  (L  -  (  -  sF- a7)+  = 


xpanding  this  in  the  usual  Taylor  series. 


31 

I  (  L' )  + 

"  3L 


0  .  1  1 

L  )  — —  + 


as 

1 

ax 


*2* 

ax 


,  32F 
1  0 


+  2  3L2 


aS.\2 

(L')|— — )  +  F ,  ( L ' .  S'.  o'.  X.  r).  t) 

o\  /  1 


3F, 

+  —  <  L'.  S'.  o'.  X.  77.  r) 


as,  aF,  as,  3F,  as, 

3X  3£  dr]  do  dr 


+  F.  = 


3F*  as,  3F*  as, 

F  (L  +  I  ,  (L  .  S  -  O  .  -.  77.  t)  + -  -  + -  -  +  F  . 

0  1  377  as  3r  do  2 
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This  can  be  separated  into  orders: 

0:  F*  =  F0(L')  =  p2/2L'2 

3Fft  a2 Fn 

-  <L'>  -  -  ,2/L'3  —  -  3/i2/L'4 

3S, 

1.  -  (m2/ L'3) -  +  F,  =  F, 

3X 

Define  F l  =  FJs  +  FJp;  ^lat  *s'  ^ i  s  ls  ^le  secular  and  F  the  periodic  part.  Then, 


* 

F,  =  F, 

1  1  S 

as,  l'3 
ax”  "  ~ijT  F,p 


The  solution  is  exactly  that  of  Brouwer  (Reference  2);  we  shall  return  to  it. 


M2 

as2 

_  _  4. 

3  H2 

/as. 

lV  + 

3F 

1  s 

as, 

3F 

4-  , 

1  S 

as, 

3F 

1  S 
-i 

3A 

2:  *  L77 

3X 

2  L'4 

\  ax 

/  + 

3L 

ax 

3* 

3  77 

3a 

a  t 

+  L( 

V  3SA 

as, 

n2 

a2s, 

as, 

n2 

a2s, 

as, 

. 

x 

C  — 

3L  \ 

L3  d\J 

ax  + 

L'3 

3£3X 

—  _j_ 

dr) 

L'3 

dodX 

3r 

1 

- 

3F 

as 

3F 

as 

1  s 

1 

1  S 

-1- 

dr) 

aT 

dr 

do 

F2  *  is  obtained  by  averaging  this  equation  over  X  to  remove  short-period  terms.  It  consists  of 
the  secular  part  of  the  summed  contributions  of  the  zonal  harmonics,  omitting  the  second, 
plus  the  solar  and  lunar  contributions,  plus  the  second-order  contribution  from  the  second  zonal 
harmonic,  which  is  the  same  as  given  in  Reference  2.  translated  into  the  canonical  variables. 
F 2  is  the  instantaneous  contribution  of  the  perturbing  potential,  except  for  the  second  harmonic. 


Dirk  Brouwer.  A.  J.  64.  378  (1959) 
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as. 


Since 


as  as 

—  =  L'  +  -  (L'.  o'.  X.  r?.  r)  +  - 

ax  ax  ax 


as, 


=  L'  + 


a2 s,  as,  a2 s,  as,  a2 s,  as,  as2 


- (L.  £.  o.  X.  rj,  r) - - -  - 

ax  3Lax  ax  a^ax  a n 


aoax  a r  ax 


we  can  solve  the  second-order  part  of  the  Hamiltonian  above: 


as. 

a2s, 

as, 

a2s, 

as, 

a2  s, 

as, 

5  /as,'' 

ax 

a  Lax 

ax- 

a^ax 

3r? 

3o3X 

a  r 

nAaTy 

L3 

L 

F,*  + 

ai- 

1  s 

as. 

ai- 

1  s 

1 

as,  +  ai^ 

as, 

+  - 

M' 

V2 

3L 

ax 

a? 

3rj  3o 

a  t 

_  ^i 

,  3S, 

9F,S 

as, ' 

3r? 

a? 

a  t 

do 

) 

Substitute  in  the  equation  for  L: 

as,  / 


L  =  L'  + 


3  as,  l3  a f 

1  1  s 


L3 


ax 


:L  ax  /a2  3L 


'  3F  as, 

+  — rl  ‘2  -  F2  +  T7 —  ^ 

M  \  3£  3r? 


ai 

as, 

ai  • 

as, 

a  i- 

as, 

1  S 

1  s 

1  S 

r)T 

a  r] 

a« 

Sr 

3o 

Sow  F  can  be  readilv  written  |  from  Brouwer  (Reference  2.  I  quation  15) 


\  ariablcs: 


I  .  (p4  R2  J,  4  F6  f2  )(-  I  +  3  p2  ) 

1  it >u i  Reference  2.  liquation  15: 


S,  K>  p‘ 


R2  J  2 / <  2 L  f)3  )  J  ”  I-  1  +  p2  >l\, 


+  Ml  +  p ),  4  L  f )  ( ( o2  -  t2  )  S  +  2o  r  ( '  )  1 

^  J  up  u  r  \ 


'  *-*  ■-*  s-’  v' s  *  ».•  *--  **• 


/* 


c.'.c  J 


NSWC  TR  84-25 


where 

Ug=  f-C  +  esinf 

Sup  =  sin  2u  +  e  sin  (u  +  c5 )  +  (e/3)  sin  (3u  -  gj) 

Cup  =  cos  2u  +  e  cos  (u  +  to)  +  (e/3)  cos  (3u  -  to) 

Appendix  D  contains  the  tedious  details  of  the  calculation  of  the  partial  derivatives  of 
S,  to  determine  the  van  Zeipel  transformation. 


as. 

as, 

aT 

x'  =  x  + 

al¬ 

as, 

as, 

dr) 

r?'  =  r?  + 

as, 

as, 

a  t 

T  =  T  + 

3a 

When  and  only  when  the  first  entry  is  made  to  the  system,  it  is  necessary  to  compute 
L’  (the  initial  value)  to  the  second  order,  in  order  that  nQ  =  n2/L '3  does  not  contain  a  second- 
order  error  (Reference  3).  For  this  case,  evaluation  of  the  derivatives  of  F  in  the  equation 
above  leads  to 


L 


=  L'  + 


as,  /  ,  as, 

aT  v  ‘  EST*  y=<s"  ’  Vl  +  1  +  2p  ’ 

.  (  as,  as,\ 

-  h,  i  +  a  +  2p-  — -  i;  — j 


5 p2 ) 


+ 


7, 


2p 


asA 

a  -  1 

dr  ) 


Note  that  this  equation  is  required  only  for  the  initial  transformation  from  osculating  to  mean 
variables  and  that  in  it  the  derivatives  of  S  are  to  be  evaluated  at  the  point  (L.  £.  a.  X.  rj.  t). 
Any  other  of  the  transformations  from  mean  to  osculating  variables  and  back  again  involves 
only  the  first  derivatives  of  S,  .  and  these  may  be  evaluated  indifferently  in  terms  of  mean  or 
osculating  variables,  since  the  difference  is  second-order  short-period  and  hence  to  be  neglected. 


K  H 


Lyddane  and  J.  Cohen. 


A.  J.  67. 


176  (1962) 
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Details  of  the  calculation  of  the  second-order  contribution  of  L'  are  found  in  Appendix.  1 


THE  HAMILTONIAN  IN  TLRMS  OF  MEAN  VARIABLES 


From  this  point  on.  all  quantities  will  be  in  terms  of  mean  variables  and  the  primes  on  the 
ariables  will  be  omitted  tor  simplicity  of  writing.  1  he  next  task  is  to  average  the  Hamiltonian 
>\ ci  the  short  period  and  to  put  the  results  into  a  form  suitable  for  calculation. 


'ONAL  HARMONICS 

We  set  aside  lor  later  consideration  the  secular  second-order  contribution  of  the  second 
■onal  harmonic  and  write  1  ,  plus  the  zonal-harmonic  portion  of  F,  as 


N  _ 

=  -  (ii,  a)  V  .1  iR  a)"  (a.r)"*1  l>  (cos  (-)) 

(Ins  comes  from  the  standard  expansion  of  an  arbitrary  axially  -s\  mmetric  gravitational  field 
nto  spherical  harmonics,  which  is  valid  as  long  as  r  >  R  .  The  bar  indicates  the  average  over 
'  i removal  of  short-period  terms),  and  (-)  is  the  colatitude  of  the  satellite,  that  is.  the  angle 
net  ween  the  unit  vector  to  the  satellite  r  and  the  unit  vector  alone  the  /-axis  r  : 

-  r 

*.  *. 

.os  (-)  -  r  •  iy 

I  he  sum  is  terminated  with  the  highest  harmonic  to  be  considered. 

Is  addition  theorem  lor  the  spherical  harmonics. 


—  (  n  m  ) 

I’  I  Os  Hi  2^,  1  -  n .  t  -  I’”1  1  COs  II  ,  )  I’"'  I  COS  d;)  COS  111  ( j  -  ig~2  ). 

"  (  ll  +  ll  l  I  1 


where  .os  (-)  ~ 
uiser'me  the  re 


.Os  II  i  COs  I) . 

presentations  ot 


+  sill  0,  sill  0 
s 

and  r  m  the 

/ 


COS  (i£,  -  I. 

A,  A  A 

L  i  1.,  Ii  trame. 


may  be  conveniently  applied  by 
namely 


v  i 


U 


S 
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(0119)  C 

(0120)  SUBROUTINE  ENTER(KEY) 

(0121)  C 

(0122)  C  KEY=  1 :  TRANSFORM  FROM  COORDINATES  TO  OSCULATING  ELEMENTS 

(0123)  C  KEY=2:  TRANSFORM  FROM  ELLIPTIC  ELEMENTS  TO  OSCULATING  ELEMENTS 

(0124)  C 

(0125)  DOUBLE  PRECISION  X ,Y , Z ,VX ,VY ,VZ , XI , ETA , SIGMA ,TAU ,L, LAMBDA ,ZETA, 

(0126)  +  L2Z,RH0,C1 ,C2, ECLON, U.COSU.SINU, ECOSF, ESINF,GAM2P,R,Z9,Z8,IC, IS, 

(0127)  +  HM,EC ,ES , A ,ECC , INC ,ELL,G ,H ,DRAD,MU 

(0128)  COMMON  /COREL/X, Y , Z , VX, VY , VZ , A ,ECC  ,INC ,ELL,G ,H/ INTEL/ETA , XI , 

(0129)  +  TAU , SIGMA , LAMBDA , L , ZET  A , L2Z , RHO , C 1 , C  2 , ECLON , U , COSU , SINU , ECOSF , 

(0130)  +  ESINF ,EC ,ES ,GAM2P/C0EF/DRAD ,MU 

(013D  C 

(0132)  IF (KEY .EQ. 1 )  GOTO  100 

(0133)  C 

(0134)  L=DSQRT (MU*A ) 

(0135)  Z8= (G+H)/DRAD 

(0136)  LAMBDA=ELL/DRAD+Z8 

(0137)  ZETA=DSQRT( 1-ECC*ECC ) 

(0138)  Z9=DSQRT(2*L*( 1-ZETA) ) 

(0139)  XI=Z9*DCOS(Z8) 

(0140)  ETA=-Z9»DSIN(Z8) 

(0141 )  Z9=2*DSIN ( INC/DRAD/2 ) *DSQRT (L*ZETA) 

(0142)  SIGMA=Z9*DC0S (H/DRAD ) 

(0143)  TAU=-Z9*DSIN(H/DRAD) 

(0144)  CALL  INTERM 

(0145)  GOTO  200 

(0146)  C 

(0147)  100  IS=-(Y*VZ-Z*VY) 

(0148)  IC=-(Z»VX-X»VZ) 

(0149)  RHO=X*VY-Y*VX 

(0150)  HMrDSQRT  ( IS**2>IC*I,2+RH0**2) 

(015D  IS  =  IS/HM 

(0152)  IC=IC/HM 

(0153)  RHO=RHO/HM 

(0154)  R  =  DSQRT(X,*2+Y**2'fZ**2) 

(0155)  Z9=(X»IS+Y»IC)/(URH0)-Z 

(0156)  U=DATAN2(Y-IC«Z9,X-IS*Z9) 

(0157)  C0SU=DC0S(U) 

(0158)  SINU=DSIN(U) 

(0159)  ECOSF =HM*HM/ (MU*R)- 1 

(0160)  ESINF= (X*VX+Y*VY+Z*VZ )*HM/ (MU*R ) 

(0161)  ZETA  =  DSQRT ( 1 -EC0SF»»2-ESINF»»  2 ) 

(0162)  EC=ECOSF*COSU+ESINF*SINU 

(0163)  ES=ECOSF*SINU-ESINF*COSU 

(0164)  Z9=ESINF/( 1+ZETA) 

(0165)  ECLON =DATAN2 ( SINU+ES-EC* Z9  »C0SU+EC+ES*Z9) 

(0166)  LAMBDA=ECLON-EC*DSIN (ECLON )  +ES*DCOS (ECLON ) 

(0167)  L=HM/ZETA 

(0168)  C1=DSQRT( ( 1+ZETA )/(2*L) ) 

(0169)  XI=EC/C1 

(0170)  ETA=-ES/C 1 

(0171)  L2Z=2*L*ZETA 
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(0107) 

KS=KS+ 1 

(0108) 

IF(MOD(KS,48).NE.3)  GOTO  200 

(0109) 

CALL  ELTRAN(KR ,2) 

(0110) 

CALL  OSTCOR 

(0111) 

C 

WRITE(1,1012)KS 

(0112) 

c 

WRITE ( 1 , 1010)A ,ECC ,INC ,ELL ,G ,H 

(0113) 

c 

WRITE ( 1,1010)(0(I), 1=1,6) 

(0114) 

c 

WRITE ( 1,1010)(VAR(I,KR),I=1,6) 

(0115) 

WRITE ( 1 , 1010)X,Y,Z,VX,VY,VZ 

(0116) 

IF(KS.LT. 1445)  GOTO  200 

(0117) 

CALL  EXIT 

(0118) 

END 
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(0107)  KS=KS+1 

(0108)  IF (MOD (KS, 48) .NE.3)  GOTO  200 

(0109)  CALL  ELTRAN(KR ,2) 

(0110)  CALL  OSTCOR 

(0111)  C  WRITE(1 ,1012)KS 

(0112)  C  WRITE ( 1 , 1010)A ,ECC , INC ,ELL,G ,H 

(0113)  C  WRITE(1,1010)(0(I),I=1,6) 

(0114)  C  WRITE ( 1,1010)(VAR(I,KR),I=1,6) 

(0115)  WRITE ( 1,1010)X,Y,Z,VX,VY,VZ 

(0116)  IF(KS.LT. 1445)  GOTO  200 

(0117)  CALL  EXIT 

(0118)  END 
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(0054) 

C 

(0055) 

C 

CYCLE  ON  STARTING  SOLUTION 

(0056) 

C 

(0057) 

120 

Z8=Z9 

(0058) 

DO  130  KR=2,KMAX 

(0059) 

Z9=0 

(0060) 

KR1=KR-1 

(0061) 

IF(KCYC.EQ.O)CALL  READ(KR) 

(0062) 

CALL  EXTRAP (KR1,KR) 

(0063) 

DO  130  N= 1 ,NEQ 

(0064) 

S=VEL(N,KR)-DEL( 1 ,N,KR) 

(0065) 

Z9=Z9+S«S 

(0066) 

DO  130  KC=1 fKR 

(0067) 

DEL(KC,N,KR)=DEL(KC,N,KR)+S 

(0068) 

130 

CONTINUE 

(0069) 

DO  150  K=2,KMAX 

(0070) 

KR=KMAX+ 1-K 

(0071) 

DO  150  N= 1 ,NEQ 

(0072) 

DO  150  KC=2,KMAX 

(0073) 

KC 1 =KMAX+2-KC 

(0074) 

150 

DEL(KC1 ,N,KR)=DEL(KC1 ,N,KR+1 )-DEL(KC1+1 ,N,KR+1 ) 

(0075) 

KCYC=KCYC+1 

(0076) 

Z9=SQRT(Z9) 

(0077) 

IF  (Z9.LT.Z8)GOTO  120 

(0078) 

IF(KCYC.LT.6)G0T0  120 

(0079) 

1010 

FORMAT (1HX,6D18.10) 

(0080) 

1012 

FORMAT (315) 

(0081) 

C 

DO  1020  K9= 1 ,KMAX 

(0082) 

K9  =  2 

(0083) 

CALL  ELTRAN(K9,2) 

(0084) 

CALL  OSTCOR 

(0085) 

C 

WRITE  (1, 1010)A,ECC,INC,ELL,G,H 

(0086) 

c 

WRITE ( 1,1010)(0(K),Kr1,6) 

(0087) 

c 

WRITE ( 1,1010)(VAR(K,K9)»K=1,6) 

(0088) 

WRITE ( 1 ,1010)X,Y,Z,VX,VY,VZ 

(0089) 

1020 

CONTINUE 

(0090) 

C 

(0091 ) 

c 

END  OF  CYCLE:  CONTINUh  THE  SOLUTION 

(0092) 

c 

(0093) 

KS=KMAX+ 1 

(0094) 

200 

KR1=MOD(KS-2,17)+1 

(0095) 

KR=MOD  (KS- 1 , 1 7 )  ■*■  1 

(0096) 

CALL  READ(KR) 

(0097) 

CALL  EXTRAP(KR1 ,KR) 

(0098) 

DO  210  N= 1 ,NEQ 

(0099) 

210 

VAR (N,KR )=VAR(N ,KR ) +ADSUM* ( VEL(N ,KR )-DEL( 1 ,N,KR) ) 

(0100) 

CALL  DERIV(KR ,  1 ) 

(0101  ) 

DO  250  N=1,NEQ 

(0102) 

S=VEL(N,KR)-DEL( 1 ,N,KR) 

(0103) 

DO  220  KC= 1 ,KMAX 

(0104) 

DEL(KC,N,KR)rDEL(KC,N,KR)>S 

(0105) 

220 

CONTINUE 

(0106) 

250 

CONTINUE 
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(0001)  DOUBLE  PRECISION  VAR( 6 , 1 7) , VEL(6 , 17) , ADAM( 17) ,DEL( 18, 6 , 17) , 

(0002)  +  STEP ,S,TIM,TIMO,ADSUM,DRAD,MU , RE,MM,MS,GAM2N ,X ,Y , Z , VX , VY ,VZ, 

(0003)  ♦  A,ECC,INC,ELL,G,H,HAR(20),0(6) 

(0004)  DIMENSION  NTERM( 3) 

( 0005 )  COMMON  // V AR , VEL , DEL , NEQ , KMAX/COEF/DRAD , MU , G AM2N , RE , MM ,MS , NTERM , 

(0006)  ♦  HAR , ADAM/COREL/X ,Y,Z,VX,VY,VZ,A, ECC , INC , ELL ,  G ,  l/,  0 

(0007)  DATA  ADAM ( 1 ) , ADAM ( 2 ) , ADAM ( 3 ) » ADAM ( 4 ) , ADAM ( 5 ) » ADAM ( 6 ) , ADAM ( 7 ) , 

(0008)  +  ADAM ( 8 ) , ADAM ( 9 ) , ADAM (10) , ADAM( 1 1 ) ,ADAM( 12) ,ADAM( 13) ,ADAM( 14) , 

(0009)  +  ADAM( 15) ,ADAM( 16) ,ADAM( 17)/ 

(0010)  ♦  -0.10000000000000D  01 , 0. 50000000000000D  00,0.83333333333333D-01, 

(0011 )  +  0.41666666666667D-01 .0.26388888888889D-01 ,0. 18750000000000D-01 , 

(0012)  +  0.1 42691 79894 180D-01, 0.11 367394179894D-01, 0. 93565365961 199D-02, 

(0013)  ♦  0. 789255401 23457D-02,0.67858499846347D-02,0.59240564123377D-02, 

(0014)  +  0. 52366932579503D- 02, 0.4677 4 984070423D- 02, 0.421 49522390055D-02, 

(0015)  ♦  0. 382689955321 20D-02, 0. 349734984 53500D- 02/ 

(0016)  DATA  HAR (2) ,HAR(3) ,HAR(4) ,HAR(5) ,HAR(6) ,HAR(7) ,HAR(8) ,HAR(9) , 

(0017)  +  HAR ( 10) ,HAR( 11 ) ,HAR( 12) »HAR( 1 3) *HAR( 14) fHAR ( 15) *HAR( 1 6) ,HAR(  17) , 

(0018)  +  HAR( 18) ,HAR( 19) ,HAR(20)/ 

(0019)  +  0.108263  4  D-02,  -0.2536  D-05, 

(0020)  +  -0.1664  D-05,  -0.2195  D-06, 

(0021)  +  0.6355  D-06,  -0.3720  D-06, 

(0022)  +  -0.3508  D-06,  -0.8733  D-07, 

(0023)  ♦  -0.5730  D-07,  0.1686  D-06, 

(0024)  +  -0.3809  D-06,  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/ 

(0025)  DATA  MU, RE, MM,MS,DRAD/. 398601  178977  8  D06,. 637814  5  D04,. 122999 

(0026)  +  7171  D-01,. 332945  561925  44  D06, .572957  795130  823  D02/ 

( 0027 )  DATA  NEQ , KMAX , STEP , TIMO , NTERM ( 1 ) , NTERM ( 2 ) , NTERM ( 3 ) / 5 , 9 , . 2 1 6D5 , 

(0028)  ♦  0. DO, 12,5,2/ 

(0029)  DATA  INIT,X,Y,Z,VX,VY,VZ,A,ECC,INC,ELL,G,H/2,0,0,0,0,0,0, 

(0030)  +  0.748503  712201  72  D04,  0.8255  D-02,  0.634300  470727  88  D02, 

(0031)  +  0.1033005  D03,  0.19952  D03,  0.1249632  D03/ 

(0032)  C 

(0033)  C  INITIALIZE 

(0034)  C 

(0035)  GAM2N= 12*HAR(2)*(MU*RE)**2 

(0036)  KCYC=0 

(0037)  Z9=1.E20 

(0038)  DO  10  1=1,17 

(0039)  DO  10  J= 1 ,6 

(0040)  DO  10  K=1,18 

(0041)  10  DEL(K,J,I)=0 

(0042)  ADSUMrO 

(0043)  DO  20  1=1, KMAX 

(0044)  ADAM ( I ) =-STEP*ADAM( I ) 

(0045)  20  ADSUM= ADSUM+ ADAM ( I ) 

(0046)  CALL  ENTER (INIT) 

(0047)  CALL  ELTRAN( 1,1) 

(0048)  CALL  READ(1) 

(0049)  CALL  DERIV (1,2) 

(0050)  DO  30  1=1,17 

(0051)  30  VAR (6 ,1) =VAR( 6 , 1 ) 

(0052)  DO  110  N= 1 ,NEQ 

(0053)  HO  DEL  (1,N,1)=VEL(N,1) 
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where  the  only  zero-order  term  is  ju2  / L3 .  a  constant.  There  are.  thus,  five  simultaneous  first- 
order  differential  equations  to  be  integrated  step  by  step  (L  being  a  constant  of  the  motion), 
after  the  partial  derivatives  of  F,  have  been  evaluated  at  each  step. 

The  program  (in  the  subroutine  DER1V)  defines  the  array 


DF(I) 


3F, 

3F, 

3  F , 

3  F , 

3F 

3a 

'  dp 

dy 

1 

1 

3a 

whose  elements  are  defined  piecemeal  in  the  previous  section  since  3/3f  =  (-  1/f2  )(3/3x).  The 
right-hand  sides  of  the  equations  of  motion  are  contained  in  the  array 


l)F<  1 ) 


3F, 

3f, 

3F , 

3F, 

3F 

,  dr] 

'  3£ 

3  r 

3a 

'  3L 

The  connection  between  these  two  is  the  Jacobian  matrix 


J(q,  p.  7-  f.  a) 

J(p.  £.  r,  o.  L) 

1'he  components  of  the  matrix  are  evaluated  in  Appendix  (1.  The  program  equates  the  matrix 
to  the  array  C(l.  J).  where  for  example  0(3.  4)  =  dyjdo.  Then  clearly. 

1)1(1)  1)1  (J  )  0(J.  I) 

j 

is  the  set  ot  time  derivatives  of  the  variables,  and  this  is  the  output  of  DFRIV. 

I  he  equations  of  motion  are  integrated  by  a  standard  Adams  numerical  integration  routine. 
Although  provision  is  made  for  an  order  as  high  as  17.  experimentation  indicated  that  an  order 
of  l)  gave  the  best  compromise  between  roundoff  error  and  truncation  error,  so.  the  program 
is  currently  set  to  l)  (KMAX  Value).  The  starting  procedure  is  found  in  the  main  program 
(lines  57-78).  At  least  (i  cycles  are  performed,  and  more  if  the  rms  change  in  the  derivatives 
per  cycle  continues  to  decline.  After  the  starting  values  are  stabilized,  the  running  procedure 
applies  a  predictor-corrector  pair  once  per  time  step.  The  extrapolation  forward  of  the  back¬ 
ward  difference  table  and  the  application  of  the  predictor  are  to  be  found  in  the  subroutine 
I  XIKAl’;  the  corrector  is  applied  in  the  main  program,  lines  l>4-108. 
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d_ 

dP 


AF. 


C 

m 


Qm 

n 


n  +  1 


THE  SOLAR  ATTRACTION 


-mS 


m  -  l 


Qm  +  I 


-n+m 

x 


(A 


m 

n  +  I 


n 


The  contribution  of  the  sun  is  exactly  the  same  as  the  moon’s  except  that  n 
are  substituted  for  pM  and  rM  and  x's,  y'g,  z^  are  substituted  for  x'M .  y'M  ,  z'M . 


THE  EQUATIONS  OF  MOTION  IN  THE  MEAN  VARIABLES 


The  Hamiltonian  in  the  mean  variables  may  be  written 
F  =  m2/2L2  +  F, 

F,  =  AF,  +  AF,  +  AFm  +  AF 

1  Z  2  s  M  S 


s  and  r< 


Since  the  variables  are  canonical,  the  Hamiltonian  equations  of  motion  are 
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THE  LUNAR  ATTRACTION 

The  moon’s  contribution  to  the  Hamiltonian  can  also  be  expanded  into  a  spherical  harmonic 
series  and  averaged  over  the  short  period: 


P  (cos  0) 


Here  cos  0  =  r.rM.  i.e..  0  is  the  angle  between  the  radius  vector  of  the  satellite  and  the  vector 
from  the  origin  to  the  moon.  The  distance  to  the  moon  is  rM  ,  and  the  components  of  the 
unit  vector  rM  are  x'M .  y'M ,  z'M .  The  expansion  is  valid  as  long  as  r  <  rM  and  the  series  is 
terminated  when  the  terms  become  negligibly  small. 


The  development  of  this  contribution  is  exactly  like  that  for  the  zonal  harmonics,  with 
two  exceptions:  a.  (3.  and  7  contain  x'M  .  y'M ,  z'M  in  place  of  x'z  y'z .  z'  and  the  coefficients 
in  x.  instead  of  Bm  (x),  are  defined  by 


-  )  cos  mf  =  em  A  m  (x) 
a  /  n  z 


The  upshot  is 


“m  /  a  \n 

iF»  s„,0i  cm «*.  mz— )  i>:  o :  it)  /r,. 

rM  m  °  "V  rM  / 


flic  sum  over  m  being  from  max  (m.  2)  by  steps  of  2  to  N. 


I  lie  recursion  relationships  for  I)"1.  0^"  (7),  and  C  (a.  j3 )  are  exactly  the  same  as  for 
the  zonal  case,  although  of  course  the  arguments  are  different.  The  recursions  and  the  derivatives 
of  the  recursive  functions  are  given  in  Appendix  F.  The  derivatives  of  AF  can  be  written  as 
follows: 


In  AF  .  obtain: 


AF„  by  replacing  C 
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THE  SECULAR  SECOND-ORDER  CONTRIBUTION  OF  THE  SECOND  HARMONIC 

Brouwer’s  Equation  29  (Reference  2)  can  readily  be  put  in  terms  of  the  variables  being 
used  here: 


A  F2s  =  (p2/L3)  y'2  2  L(f/24)  }  S  f 2  (1  -  18  p2/ 5  +  p4) 

+  4  f  (1-6  p2  +  9p4 )  -5  (1-2  p2  -7  p4)  -  2(1-15  p2)(a2  -  /32)j 
It  follows  that,  if  M  =  (p2/L3)  yn2  L  f/24  and  7  =  p, 

7-  AF  =  M(  -  4a)(l- 1 5  p2) 

3a  2S 

0 

—  A  F2$  =  M(4  jJXI-15  P2) 

^  AF2s  =  M(4  p)  j  f2  (-9  +  5 p2)  +  12  f  (-1  +  3  p2) 

+  5  (1  +  7  p2)  +  15  (a2  -  02)| 

0 

—  AF,s  =  M(-  J2)  j  10  f  (1-18  p2  /  5  +  p4)  +  4  (1-6  p2  +  9  p4 )  ( 

3x  2S  '  > 

+  7  f  •  AF 

J  2  s 

—  AF,.  =  -  (5/a)  AF,. 
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Then,  a  =  ec  ^  +  e$ 

0  =  es  i/'i  -  ec  i^2 

and  Cm  (a,  0)  =  Re  [(a  +  i  0)m  |  =  em  (1  -  72)m/2  cos  m  (gj  -  <£2 ) 
Introduce  tor  convenience  one  more  set  of  polynomials  by 


I)'” 

n 


( n-  m )  ! 

- —  Qm 

(n+m)  !  n 


(0) 


and  introduce  the  new  variable  \  =  1/f. 


1  lien. 


AF  = 


<-  -  6,no»  ('n,  «*■  <M/a)  (Re/a)n  D"’  Q”  (7)  B„m 


(X) 


The  sum  over  11  starts  at  2  or  111.  whichever  is  greater,  and  goes  in  steps  of  2  to  N. 

This  formulation  lias  the  advantages  that  the  Hamiltonian  depends  only  on  five  intermediate 
variables  [a.  0.  7.  f  (or  x),  and  a|  through  three  functions  <Q™.  B™ .  and  C'  )  and  one  set 
of  constants  ( l)"’  ).  all  of  which  can  be  calculated  recursively:  and  the  partial  derivatives  are 
reasonably  simple. 

I  he  recursion  relationships  and  the  derivatives  of  the  functions  arc  developed  in  Appendix  T. 

I  he  derivatives  of  AT  can.  with  the  aid  of  those  results,  be  written  as  the  following  schema: 

/ 

In  AT  .  obtain. 

/ 


1) 

—  A  1  In  replacing  ('  with  mC 
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A  ^  A  S 

rz:  sin  02  cos  «p2  =  ^,  =  L,  *  rz  =  x'  +  —  (7  +  z[) 

A  lc 

sin  02  sin  <p2  =  2  =  L2  *f  =  y'  +  —  (7  +  z') 

cos  02  =  7  =  h*  r  =  -  i  x'  -  1  y '  +  p  z' 

*  1  z  sz  c  z^z 

where  the  unit  vector  r 

(There  is  an  obvious  simplification  if  the  actual  values  for  xz,  y'.  z'  are  substituted  in  the 
foregoing,  but  this  form  is  used  for  parallelHsm  with  the  lunar  and  solar  contributions  to  the 
Hamiltonian.) 

Pn  (cos  0)  =  £o  (2  -  8m0)  K  <°)  ?n  (7)  cos  m  (u  -  *>2) 

If  Q™  (7)  =  dm/d7m  Pn  (7),  then  P™  (7)  =  (1  -  y2)ml2  Q™  (7)  and 

-1^  f  j-j—  j-j-j  \  f 

(a/r)n  + 1  P  (cos  0)  =L  (2  -  5  )  - - 7  Qm  (0)  Qm  (7)  (a/r)n+l  cos  mf- 

m  mo  (n+m)  !  "  "  ' 

(T  -  y2)m  / 2  cos  m  (cb  -  ip2 ) 

We  define  B™  by  (a/r)n+1  cos  mf  =  em  B™  (knowing  that  the  Cayley  coefficients  are 
proportional  to  em )  and  consider 


[e  v7 1-72  exp  i  (a>  -  <p2)lm  =  l(ec  +  i  %)  (>Pi  ~  i  (MJm 
=  (ec  (//]  +  es  \p2  +  i  (es  1 //,  -  ec  ^2)]m 

Define  polynomials  Cm  and  Sm  by 

(a  +  i  /3)m  =  Cm  (a,  0)  +  i  Sm  (a,  0) 
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(0172)  C2=DSQRT(( 1+RH0)/L2Z) 

(0173)  SIGMA=IC/C2 

(0174)  TAUsI S/C2 

(0175)  C 

(0176)  200  CALL  XFR(0,3) 

(0177)  RETURN 

(0178)  END 
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(0179)  C 
(0180) 
(0181)  C 
(0182)  C 
(0183)  C 
(018*1) 
(0185) 
(0186) 
(0187) 
(0188) 
(0189) 
(0190)  C 
(0191 ) 
(0192) 
(0193)  C 
(019*0 
(0195) 
(0196) 
(0197) 
(0198) 
(0199) 
(0200) 
(0201  ) 
(0202) 
(0203) 
(020*0  C 
(0205) 
(0206) 
(0207) 
(0208) 
(0209) 
(0210) 
(0211  ) 
(0212) 
(0213) 
(021*0 


SUBROUTINE  OSTCOR 

TRANSFORM  FROM  OSCULATING  ELEMENTS  TO  COORDINATES/ELLIPTIC  ELEMENTS 

DOUBLE  PRECISION  X ,Y , Z ,VX ,VY ,VZ ,XI , ETA , SIGMA ,TAU ,L .LAMBDA ,ZETA, 

+  L2Z,RH0,C1 ,C2,ECLON,U,COSU,SINU,ECOSFtESINF,R,Z9,EC,ES, 

♦  A ,ECC , INC , ELL,G ,H ,DRAD ,MU 

COMMON  /COREL/X ,Y ,Z ,VX , VY , VZ , A ,ECC , INC , ELL, G ,H/INTEL/ETA ,XI , 

+  TAU , SIGMA , LAMBDA ,L,ZETA ,L2Z ,RHO,C 1 ,C2,ECL0N ,U ,COSU ,SINU ,ECOSF , 

+  ESINF,EC,ES/COEF/DRAD,MU 

CALL  XFR(0, 1 ) 

CALL  INTERM 

A=L*L/MU 

ECC=DSQRT( 1-ZETA*ZETA) 

Z9=DSQRT( 1-RH0*RH0) 

INC=DRAD»DATAN2(Z9,RH0) 

IF(SIGMA.NE.O.OR.TAU.NE.O)  H = DR AD*DATAN2 ( -T AU , SIGMA ) 

G=DRAD*DATAN2(-ETA ,XI )-H 

ELL=DRAD*LAMBDA-G-H 

ELL=DMOD(ELL,3.6D2) 

IF ( ELL . LT . 0 . DO ) ELL  =  ELL+  3 • 6D2 
IF(ELL.GT. 1 . 8D2)ELL=ELL-3.6D2 

R= ( 1-EC*DC0S(ECL0N )-ES*DSIN(ECLON) )*L*L/MU 
X=R* (COSU« ( 1 -TAU*TAU/L2Z)-SINU*SIGMA*TAU/L2Z ) 

Y=R« (-C0SU»SIGMA»TAU/L2Z+SINU» ( 1-SIGMA»SIGMA/L2Z ) ) 
Z9=DSQRT((URH0)/L2Z) 

Z=K*Z9* (SINU»SIGMA+COSU«TAU ) 

VX=((ECL0N-LAMBDA)«X+ZETA«(-C2*SIGMA"Z-FH0»Y))«L/(R«R) 

VY=((ECL0N-LAMBDA)«Y+ZETA*(RH0»X+C2«TAU»Z))»L/(R»R) 

VZ=((ECL0N-LAMBDA)«Z^ZETA»(-C2*TAU»Y+C2»SIGMA»X))«L/(R»R) 

RETURN 

END 
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(0215)  C 
(0216) 
(0217)  C 
(0218)  C 
(0219)  C 
(0220)  C 
(0221 ) 
(0222) 
(0223) 
(0224) 
(0225) 
(0226) 
(0227) 
(0228)  C 
(0229) 
(0230)  20 
(0231) 
(0232) 
(0233) 
(0234) 
(0235) 
(0236) 
(0237) 
(0238) 
(0239) 
(0240) 
(0241) 
(0242) 
(0243) 
(0244) 
(0245) 
(0246) 
(0247) 
(0248) 
(0249) 
(0250) 
(0251) 
(0252) 
(0253) 
(0254) 
(0255) 
(0256) 
(0257) 
(0258) 
(0259) 
(0260) 
(0261) 
(0262) 
(0263) 
(0264) 
(0265) 
(0266)  C 
(0267) 


SUBROUTINE  ELTRAN (KROW ,KEY ) 

KEY=1:  TRANSFORM  FROM  OSCULATING  TO  MEAN  ELEMENTS 
KEY=2:  TRANSFORM  FROM  MEAN  TO  OSCULATING  ELEMENTS 

DOUBLE  PRECISION  XI , ETA .SIGMA ,TAU ,L .LAMBDA ,ZETA, 

+  L2Z.RH0.C1 ,C2,ECLON,U,COSU,SINU,ECOSF,ESINF,GAM2P,C(6) ,D(6) , 
+  E(6),0(6),ZETA2,RH03,Z9,Z8,Z7,Z6,Z5,Z4,Z,C2U,S2U,SU3,CU3,UL, 
+  SUP ,CUP .SUPS2.SGH ,EC ,ES , VAR ( 6 , 1 7 ) ,DRAD ,MU 
COMMON  /COREL/C ,E,0// VAR/ INTEL/ETA, XI ,TAU .SIGMA .LAMBDA ,L , 

+  ZETA ,L2Z ,RHO,C 1 ,C2, ECLON ,U ,COSU ,SINU .ECOSF .ESINF ,EC ,ES .GAM2P 
+  /COEF/DRAD.MU 

CALL  XFF (KROW, KEY) 

CALL  INTERM 
ZETA2=ZETA*ZETA 
Z= 1+ECOSF 

RH03=(-U3*RH0»RH0)/3 

Z9=SIGMA*SIGMA-TAU*TAU 

Z8=2*SIGMA*TAU 

Z7=DSIN(2*U) 

Z6=DCOS(2*U) 

Z5=SINU+DSIN ( 3*U )/3 
Z4=COSU+DCOS ( 3*U ) / 3 
C2U=Z6*Z9-Z7*Z8 
S2U=Z7*Z9+Z6*Z8 
SU3=Z5*Z9+Z4*Z8 

CU3= (2*COSU-Z4)*Z9+ (-2*SINU+Z5)*Z8 
UL=U-LAMBDA+ESINF 
SUP=Z7+EC*Z5+ES* (2*COSU-Z4) 

CUP=Z6+EC*Z4+ES* (Z5-2*SINU) 

SUPS2=SUP*Z9+CUP*Z8 

SGH  = ( 1 +RHO* ( 2-5*RHO ) ) *UL- ( 3+5«RHO )«SUPS2/ ( 2*L2Z ) 

Z9=Z* ( 1+Z)* (RH03+ ( 1 +RHO)/L2Z*C2U ) 

D(5) =-GAM2P* (SGH- (ESINF* Z9+ZETA2* (RH03*ESINF+( 1+RH0)/(2*L2Z)* 
♦  (SUPS2-S2U) ) )/( 1+ZETA) ) 

Z8= (Z/ZETA)**3 

D(6) =GAM2P* (RH03* (Z8- 1 )*L2Z+ ( 1+RHO)*Z8*C2U )/2 
Z9=(RH03+Z*(URH0)*C2U/L2Z)/(2»ZETA2) 

Z8= ( 1+Z)*Z9+RH03* (Z/ZETA)**2/2 

Z7=ESINF*( 1+2*ZETA)*Z8+ ( 1+RH0)*(SUPS2-S2U)/(8*L) 

Z7=SGH+RH03+Z7/ ( 1+ZETA) 

Z6= ( 1+ZETA+ ZETA2)*Z9 
Z5= ( 1+RH0)*C 1/4 
Z9=2*L*C 1 *Z8 

D(  1  )=-GAM2P*(-XI*Z7+Z9*SINU-Z6*ETA+Z5*SU3) 

D(2)=GAM2P* (-ETA*Z7+Z9*COSU+Z6*XI-Z5*CU3) 
Z9=2*RH0*UL+SUPS2/(2*L2Z) 

D(3)=-GAM2P*(-SIGMA*Z9+(URHO)/2»(SIGMA«SUP+TAU*CUP)) 

D( 4) =GAM2P» (-TAU»Z9+ ( 1 +RH0)/2» (-TAU*SUP+SIGMA*CUP) ) 

IF(KEY.EQ. 1 )GOTO  100 


i 


1 

J 


i 
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(0268)  C 

(0269)  DO  50  1*1,6 

(0270)  50  0(I)=VAR(I,KR0W)+D(I) 

(0271)  RETURN 

(0272)  C 

(0273)  100  Z8=1+RHO«(2-5»RHO) 

(0274)  D(6)=D(6)«(1-3/(2»L)«D(6)+GAM2P»(-3»ZETA*RH03+Z8))+ 

( 0275)  ♦  GAM2P* (-Z8* (ETA*D ( 1 ) +XI*D ( 2 ) )-2»RHO» (TAU*D ( 3 )+SIGMA*D ( 4 ) ) ) 

(0276)  DO  150  1*1,6 

(0277)  150  VAR(I,KR0W)=0(I)-D(I) 

(0278)  RETURN 

(0279)  END 
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(0280)  C 
(0281) 
(0282)  C 
(0283)  C 
(0284)  C 
(0285) 
(0286) 
(0287) 
(0288) 
(0289)  C 
(0290) 
(0291) 
(0292) 
(0293) 
(0294) 
(0295) 
(0296) 
(0297) 
(0298) 
(0299)  100 
(0300) 
(0301) 
(0302) 
(0303) 
(0304) 
(0305) 
(0306) 
(0307) 
(0308) 
(0309) 
(0310) 
(0311) 


SUBROUTINE  INTERM 

CALCULATE  INTERMEDIATE  QUANTITIES  IN  ELEMENT  TRANSFORMATIONS 

DOUBLE  PRECISION  XI , ETA, SIGMA ,TAU,L, LAMBDA, ZETA.L2Z.RH0.C 1 ,g2, 
*  ECLON ,U ,C0SU ,SINU .ECOSF ,ESINF ,EC ,ES , DRAD ,MU ,GAM2N ,GAM2P , Z8 , Z9 
COMMON  /INTEL/ETA , XI ,TAU , SIGMA .LAMBDA ,L , ZETA ,L2Z , RHO ,C 1 ,C2 , 

+  ECLON ,U ,COSU ,SINU , ECOSF .ESINF , EC ,ES .GAM2P/C0EF/DRAD ,MU .GAM2N 

ZETAs 1- (XI*XI+ETA*ETA ) / ( 2«L ) 

L2Z=2*L*ZETA 
GAM2P=GAM2N/L2Z*»4 
RHO= 1-(SIGMA*SIGMA+TAU*TAU)/L2Z 
C 1=DSQRT( ( 1+ZETA)/ (2*L ) ) 

C2rDSQRT( ( URH0)/L2Z) 

EC=C1«XI 

ES=-C1*ETA 

ECLON=LAMBDA 

Z9=ECL0N 

Z8=1-EC»DcoS(Z9)-ES»DSIN(Z9) 

ECL0N=Z9- ( Z9-LAMBDA-EC*DSIN ( Z9 ) +ES»DCOS ( Z9 ) ) / Z8 
IF(DABS(ECL0N-Z9).GT.2.D-14»(DABS(Z9)+1))G0T0  100 
Z9= (ECLON-LAMBDA ) 

Z8=Z8+ZETA 

U  =  ECLON+2*DATAN2(Z9»  Z8) 

SINU=DSIN (U ) 

COSU=DCOS(U ) 

ESINF =EC*SINU-ES*COSU 

ECOSF=EC*COSU+ES*SINU 

RETURN 

END 
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(0312) 

C 

(0313) 

SUBROUTINE  XFR (KROW, KEY) 

(0314) 

C 

(0315) 

C 

TRANSFER  OSCULATING  ELEMENTS  0  TO  WORKING  ELEMENTS  W  (KEY  - 1 ) 

(0316) 

C 

OR  MEAN  ELEMENTS  VAR  TO  WORKING  ELEMENTS  (KEY=2)  OR  THE  REVERSE 

(0317) 

C 

(KEY =3  OR  4). 

(0318) 

C 

(0319) 

DOUBLE  PRECISION  C(6) ,E(6) ,0(6) ,VAR(6, 17) ,W(6) 

(0320) 

COMMON  / INTEL /W/COREL/C ,E ,0//VAR 

(0321  ) 

C 

(0322) 

GOTO  (10, 20, 30, 40), KEY 

(0323) 

10 

DO  15  1=1,6 

(0324) 

15 

W(I)=0(I) 

(0325) 

RETURN 

(0326) 

20 

DO  25  1=1,6 

(0327) 

25 

W(I ) =VAR( I ,KROW) 

(0328) 

K9=KR0W 

(0329) 

RETURN 

(0330) 

30 

DO  35  1=1,6 

(033D 

35 

0(I)=W(I) 

(0332) 

RETURN 

(0333) 

40 

DO  45  1=1,6 

(0334) 

45 

V  AR ( I , KROW ) =W ( I ) 

(0335) 

RETURN 

(0336) 

END 
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0334 
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< 

(0337) 

C 

(0338) 

SUBROUTINE  READ (KROW) 

* 

(0339) 

C 

1 

(0390) 

c 

READ  SUN  AND  MOON  COORDINATES  FROM 

DISK  FILE 

(0391) 

c 

* 

(0392) 

DOUBLE  PRECISION  TSL ,SL( 1 7 , 8) , ZZ 

• 

(0393) 

COMMON/LUNSOL/TSL ,SL 

(0399) 

READ(5 , 100)TSL, (SL(KROW,I ) ,1=5 ,7) 

. 

(0395) 

READ(5, 100) (SL(KROW, 

I), 1=1, 3) 

• 

(0396) 

READ( 5 , 1 00 ) ZZ 

(0397) 

READ ( 5 , 1 00 ) ZZ 

(0398) 

100  FORMAT (4D22. 19) 

(0399) 

DO  200  1=1,2 

(0350) 

N=9» (1-1 ) 

! 

(0351) 

SL(KROW ,N+4 )  =  DSQRT  (SL(KR0W,N+1  )  #*2+SL(  KROW  ,N  +  2 )  *,*2+SL  (KROW . 

(0352) 

+  N+3)**2) 

(0353) 

DO  150  J=  1 , 3 

(0359) 

150  SL(KROW ,N+ J ) =SL (KROW 

,N+J)/SL(KROW,N 

+  4 ) 

(0355) 

200  CONTINUE 

(0356) 

RETURN 

• 

(0357) 

END 
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(0358) 

C 

(0359) 

SUBROUTINE  QUE (M ,Q0 ,Y ,NT ,Q , J ) 

(0360) 

C 

(J361 ) 

C 

Q(JS,N)=  Q  M/N  :  Q(JL,N)=  Q  M+1/N 

(0362) 

C 

(0363) 

DOUBLE  PRECISION  Q0(2,20),Y,Q 

(0364) 

IF (M.GT. 1 )  Q0(J ,M- 1 )=0 

(0365) 

IF (M .GT .NT ) RETURN 

(0366) 

IF(M.EQ.0)G0T0  100 

(0367) 

Q0( JtM)=Q 

(0368) 

IF (M. EQ. NT ) RETURN 

(0369) 

Q0( J ,M+ 1 )= (2*M+1 )*Y*Q 

(0370) 

N=M+2 

(037D 

GOTO  200 

(0372) 

100 

Q0( J , 1 )=Y 

(0373) 

Q0( J ,2) = 1 . 5*Y*Y-. 5 

(0374) 

N  =  3 

(0375) 

150 

Q0(J,N)=(( 2*N- 1 )*Y*Q0( J ,N-1 )- (N+M-1 )*Q0( J ,N-2) )/ (N-M) 

(0376) 

N=N+ 1 

(0377) 

200 

IF(N.LE.NT)GOTO  150 

(0378) 

Q=(2*M+1 )*Q 

(0379) 

RETURN 

(0380) 

END 
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(0301) 

c 

(0382) 

SUBROUTINE  AYBEE(M,B0 ,X ,NT , B , J , IDEC ) 

(0383) 

c 

(0384) 

DOUBLE  PRECISION  B0(2,20) , B,X,X2,P0,P1 

(0385) 

X2=X*X 

(0386) 

IF(IDEC.EQ. 1 )GOTO  200 

(0387) 

c 

(0388) 

c 

B0(JS,N)=  A  M/N+2  :  B0(JL,N)=  A  M+1/N+2 

(0389) 

c 

(0390) 

P0=B 

(0391  ) 

P 1 = ( 2*M> 1 )*B/ (M+ 1 ) 

(0392) 

B=-P  1 

(0393) 

N=M+ 1 

(0394) 

IF(M.GT.2)  B0(J,M-2)=P0 

(0395) 

IF(M.GT.I)  BO ( J ,M- 1 ) =P1 

(0396) 

GOTO  150 

(0397) 

100 

B0( J,N)s ( (2*N+ 1 )*P1- (N+M)* (N-M)*P0/ (X2*N ) )/ (N+ 1 ) 

(0398) 

PO=P1 

(0399) 

P1=B0(J,N) 

(0400) 

N=N+ 1 

(0401  ) 

150 

IF (N . GT. NT ) RETURN 

(0402) 

GOTO  100 

(0403) 

C 

(0404) 

C 

B0(JS,N)=  B  M/N  :  BO(JL,N)=  B  M+1/N 

(0405) 

C 

(0406) 

200 

IF(M.EQ.O)  GOTO  210 

(0407) 

B0(J,M)=0 

(04U8) 

IF(M.GT.I)  B0(J,M-1)=0 

(0409) 

IF (M.GE .NT ) RETURN 

(0410) 

210 

B0(J,M+1 )sB 

(0411) 

IF (M .EQ.NT- 1 ) RETURN 

(0412) 

B0(J,M+2)=(M+1 )*X2*B 

(0413) 

N=M>3 

(0414) 

GOTO  250 

(0415) 

220 

B0(J,N)=(N-1 )*X2/( (N-M-1 )*(N+M-1 ))*((2"N-3)*B0(J,N-1) 

(0416) 

♦  *B0( J ,N-2) ) 

(0417) 

N=N+ 1 

(0418) 

250 

IF(N.LE.NT)GOTO  220 

(0419) 

B=B*X2/2 

(0420) 

RETURN 

(0421  ) 

END 
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B 

D 

ARGUMENT 

000007 

0382S 

0419M 

0384S 

0390 

0391 

0392M 

0410 

0412 

BO 

D 

ARGUMENT 

000004 

0382S 

0384S 

0394M 

0395M 

0397M 

0399 

0407M 

0408M 

0410M 

0412M 

0415M 

IDEC 

I 

ARGUMENT 

00001 1 

0382S 

0386 

J 

I 

ARGUMENT 

000010 

0382S 

039*1 

0395 

0397 

0399 

0407 

0408 

0410 

0412 

0415 

M 

I 

ARGUMENT 

000003 

0382S 

0391 

0393 

0394 

0395 

0397 

0406 

0407 

0408 

0409 

0410 

0411 

0412 

0413 

0415 

N 

I 

000371 

0393M 

0397 

0399 

0400M 

0401 

041  3M 

0415 

0417M 

0418 

NT 

I 

ARGUMENT 

000006 

0382S 

0401 

0409 

0411 

0418 

PO 

D 

000372 

0384S 

0390M 

0394 

0397 

0398M 

PI 

D 

000376 

03845 

0391M 

0392 

0395 

0397 

0398 

0399M 

X 

D 

ARGUMENT 

000005 

0382S 

0384S 

0385 

X2 

D 

000402 

0384S 

0385M 

0397 

0412 

0415 

0419 

100 

000106 

0397D 

0402 

150 

000172 

0396 

040 1 D 

200 

000177 

0386 

0406D 

210 

000227 

0406 

0410D 

220 

000263 

0415D 

0418 

250 

000344 

0414 

0418D 

NSVVC  TR  84-25 


TC  OR 


By  .1  cal!  to  XI  R.  mines  the  osculating  elements  into  the  working  area.  Call  IN  1 1  RM  to 
get  the  intermeiliate  quantities  appropriate  to  those  elements.  Calculate  both  elliptic  elements 
and  (  artesian  coordinates. 


IRAN  (RROVV.  KLY) 

KROW  points  to  the  step  value  in  VAR  to  which  the  mean  elements  apply 
kl  Y  -  I 


Call  XI  R  and  IN  11  RM  to  move  the  osculating  elements  into  position.  Calculate  the  short- 
period  terms  that  are  the  differences  between  mean  and  osculating  elements.  Apply  these 
and  store  the  mean  elements  obtained  in  the  proper  place  in  VAR. 

NO  IT  :  In  this  process,  only  part  of  the  second-order  short-period  correction  to  L  is  made, 
the  remainder  being  done  more  conveniently  in  DCRIV. 

kl  Y  -  2 

I  \ecute  the  reverse  process,  from  mean  elements  to  osculating  elements. 


It  KM 


Csing  the  working  elements  as  a  base,  calculate  the  intermediate  quantities  (e.g..  ZI  TA. 
I  5INI  .  1C)  (Here  also  Kepler's  equation  is  solved). 


R 


Moves  osculating  elements  or  mean  elements  into  working  area  or  reverse. 


AO 


Reads  the  file  L  I  NSOL  to  acquire  lunar  solar  coordinates  and  calcu.ates  the  unit  vector 
and  the  distance  to  each  body. 


4 


1 


A 


"  a 


B-4 


-4  4^-^- 
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GPS 

1.  The  program  starts  witli  DATA  statements  defining  the  numerical  constants  and  the  initial 
parameters  of  the  problem,  which  include  the  initial  Cartesian  coordinates  or  elliptic  elements 
for  the  satellite  in  the  last  of  the  DATA  statements  (line  29  of  the  listings). 

2.  After  initialization  and  calculation  of  the  starting  values  of  the  mean  elements  and  their 

time  derivatives,  the  program  in  line  57  begins  the  cycling  to  set  up  the  starting  solution. 

3.  When  this  is  complete,  lines  81  through  89  allow  printing  information  about  the  starting 

solution.  Note  that  all  of  these  lines  can  be  commented  out  or  removed  without  affecting 
the  operation  of  the  program.  See  also  the  note  below  in  paragraph  5. 

4.  Lines  93  through  107  then  apply  one  predictor  and  one  corrector  step  to  move  the  solution 
forward  one  time  step:  line  116  governs  the  repetitions  of  this. 

5.  Lines  108  through  115  allow  printing  information  about  the  continuing  solution,  again 

these  lines  can  be  changed  or  removed. 

NOTE:  The  form  of  the  write  statements  in  this  program,  as  well  as  that  of  the  READ 
statements  in  subroutine  READ  (lines  344-347).  might  need  to  be  altered  for  running  the 
program  on  a  different  computer.  PRIME  Fortran  IV  (in  which  the  program  is  written) 
recognizes  logical  unit  1  as  the  terminal  and  logical  unit  5  as  the  first  file  channel.  The 
proper  file  for  reading  is  opened  at  the  operating  system  level  by  OPEN  LUNSOL  1  1. 

which  opens  the  disk  file  LUNSOL  on  file  channel  1  for  reading  only. 


ENTER  (KEY) 

1.  KEY  =  1 

Transforms  Cartesian  coordinates  (position  and  velocity)  into  osculating  elements  and.  by  a 
call  to  INTERM,  prepares  the  intermediate  quantities  useful  in  further  transformations. 
Calls  XFR  to  move  the  working  set  of  osculating  elements  to  safe  storage. 

2.  KEY  =  2 

Transforms  elliptic  elements  into  osculating  elements  (calculating  the  intermediate  quantities 
in  the  process)  and.  by  calling  XFR.  moves  the  working  elements  into  safe  storage  as  the 
osculating  elements. 
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APPENDIX  B 

PROGRAM  AND  SUBROUTINES 


B-l 
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(0656)  C 
(0657) 
(0658)  C 
(0659) 
(0660) 

(0661 ) 
(0662) 
(0663) 
(0664) 
(0665) 
(0666) 
(0667) 
(0668)  100 
(0669)  110 

(0670) 

(0671  ) 
(0672) 


SUBROUTINE  EXTRAP (KROW 1 ,KROW ) 

DOUBLE  PRECISION  DEL ( 1 8,6 , 1 7) , ADAM( 1 7) , VAR( 6 , 1 7 ) , VEL (6 , 1 7 ) , 
♦  S,XX(6) ,HH(20) 

DIMENSION  NT( 3) 

COMMON  //VAR , VEL ,DEL ,NEQ .KMAX/COEF/XX ,NT ,HH , ADAM 
DO  110  N= 1 ,NEQ 
S  =  0 

DO  100  KC=  1 ,KMAX 
KC 1 =KMAX+ 1-KC 

DEL(KC 1 ,N ,KROW) =DEL (KC 1 ,N ,KR0W1 )+DEL (KC 1 + 1 ,N ,KROW ) 
S=S+ADAM(KC1 )«DEL(KC1 ,N,KROW) 

VAR (N ,KROW ) =VAR  (N ,KROW 1  )+S 
CALL  DERIV (KROW, 1  ) 

RETURN 

END 


ADAM 

D 

/COEF/ 

000153 

0659S 

0662S 

0668 

DEL 

D 

// 

001460 

0659S 

0662S 

0667M 

*  0668 

DERIV 

R 

EXTERNAL 

000000 

0670 

KC 

I 

000146 

0665M 

0666 

KC  1 

I 

000147 

0666M 

0667 

0668 

KMAX 

I 

// 

017741 

0662S 

0665 

0666 

KROW 

I 

ARGUMENT 

000004 

0657S 

0667 

0668 

0669  0670A 

KR0W1 

I 

ARGUMENT 

000003 

0657S 

0667 

0669 

N 

I 

000150 

0663M 

0667 

0668 

0669 

NEQ 

I 

// 

017740 

0662S 

0663 

S 

D 

000151 

0659S 

0664M 

0668M 

0669 

VAR 

D 

// 

000000 

0659S 

0662S 

0669M 

100 

000053 

0665 

0668D 

110 

000106 

0663 

0669D 

A -30 
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002411 

0484 

0510D 

002537 

0517 

0518 

051 

002662 

0533 

0535D 

002702 

0529 

0537D 

002752 

0542 

0543D 

002771 

0536 

0544D 

003037 

0551 

0553D 

003115 

0557 

0562D 

003175 

0561 

0569D 

003246 

0573 

0576D 

003452 

0577 

0588 

055 

003525 

0594 

0596D 

003571 

0555 

0576 

059 

003577 

0597 

0599D 

004165 

0599 

0617D 

004171 

0617 

0618 

061 

004350 

0628 

0630D 

004401 

0623 

0632D 

004547 

0641 

0643D 

004572 

0443 

0620 

063 

004637 

0645 

0648D 

001170 

0438 

0440D 

001434 

0465 

0466 

046 

002006 

0489 

0490 

049 
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QUE 

R  EXTERNAL 

000000 

0558 

0570 

0627 

0640 

R 

D 

005027 

0428S 

0622M 

0624 

0625 

0632 

0635 

J 

RE 

D  /COEF/ 

000014 

0428S 

0434S 

0531 

0624 

RHO 

D  /INTEL/ 

000040 

0428S 

0434S 

0451 

0452 

0458 

0459 

0468 

% 

0483 

0503 

0504 

0509 

0603 

0612 

0636 

RH02 

D 

005033 

0428S 

0603M 

0605 

0606 

0608 

0609 

0612 

' 

SIGMA 

D  /INTEL/ 

000014 

0428S 

0434S 

0449 

0450 

0454 

0455 

0460 

0462 

0471 

0485 

0494 

0498 

0507 

0621 

*1 

0637 

SINB 

D 

005037 

0428S 

0621M 

0627A 

0636 

0637 

SINU 

D  /INTEL/ 

000070 

0428S 

0434S 

0621 

0637 

SL 

D  /LUNSOL/ 

000004 

0428S 

0434S 

0485 

0486 

0487 

0488 

0491 

0493 

0494 

0495 

0496 

0498 

0499 

0500 

* 

0501 

0507 

0508 

0537 

0632 

0637 

SM 

D 

005043 

0428S 

0550M 

0566 

0568M 

0597 

SMI 

D 

005047 

0428S 

0566M 

0567 

0568 

0592 

TAU 

D  /INTEL/ 

000010 

0428S 

0434S 

0449 

0450 

0453 

0456 

0461 

0463 

0472 

C485 

0493 

0499 

0508 

0621 

0637 

V 

D 

005053 

0428S 

0609M 

0614 

0615 

0616 

- 

VI 

D 

005057 

0428S 

0449M 

0451 

0481 

0486 

0493 

0494 

0495 

0496 

0497 

4 

V2 

D 

005063 

04282. 

04'j0M 

0458 

0482 

0487 

0498 

0499 

; 

0500 

0501 

0502 

VAR 

D  // 

000000 

0428S 

0434S 

0647M 

VEL 

D  // 

000630 

0428S 

0653M 

0434S 

0648M 

0649M 

0650M 

0651M 

0652M 

X 

D 

005067 

0428S 

0527M 

0548 

0559A 

0571A 

0585 

XFR 

R  EXTERNAL 

000000 

0441 

XI 

D  /INTEL/ 

000004 

0428S 

0434S 

0449 

0450 

0453 

0456 

0460 

4 

0462 

0469 

0475 

0486 

0487 

0493 

0496 

4 

0498 

0500 

0505 

0511 

0514 

J 

Z5 

D 

005073 

0428S 

0608M 

0609 

0610 

0611 

. 

Z6 

D 

005077 

0428S 

0607M 

0609 

0612 

Z7 

D 

005103 

0428S 

0485M 

0486 

0487 

0488 

0492M 

0493 

0494 

0495 

0496 

0497 

0498 

0499 

0500 

0501 

0502 

0503 

0538M 

0539M 

0540 

0656M 

' i 

0609 

0614 

0625M 

0629M 

0630 

0633M 

0634M 

0635M 

0642M 

0643 

Z8 

D 

005107 

0428S 

0451M 

0452 

0455 

0456 

0458M 

0459 

0462 

0463 

0504M 

0505 

0506 

0532M 

0534M 

0535 

0540M 

0541 

0582M 

0583 

0586 

0589 

0605M 

0609 

0614 

0624M 

0625 

0629 

0632M 

0635 

0642 

Z9 

D 

003113 

0428S 

0452M 

0453 

0454 

0457 

0459M 

0460 

0461 

0464 

0468M 

0469 

0470 

0473 

0474M 

0475 

0476 

0479 

0503M 

0504 

0507 

0508 

• 

0509 

051 0M 

0511 

0512 

0513 

0514 

0515 

0516 

0531M 

0532 

0534 

0537M 

0540 

0541 

0543 

0581M 

0582 

0585 

0604M 

0609 

0610 

\ 

0611 

0612 

0614 

0636M 

0637M 

0640A 

ZETA 

D  /INTEL/ 

000030 

0428S 

0434S 

0474 

0479 

0510 

0513 

0516 

' 

05 2/ 

0605 

0609 

0612 

0614 

0622 

; 

A-28 

i 

| 

V  W  V'.'  VvV*/* 

•*-  •*.  -* 
-■  'U  -- 

i 
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ECOSF 

D 

/INTEL/ 

000074 

0428S 

0434S 

0622 

ETA 

D 

/INTEL/ 

000000 

0428S 

0434S 

0449 

0450 

0454 

0455 

0461 

0463 

0470 

0476 

0486 

0487 

0494 

0495 

0499 

0501 

0506 

0512 

0515 

FO 

D 

000546 

0428S 

0552M 

0583M 

0584M 

0585M 

0586M 

0591 

0592 

0593M 

0595 

0596M 

GAM2P 

D 

/INTEL/ 

000114 

0428S 

0434S 

0604 

GAMMA 

D 

004775 

0428S 

0483M 

0488M 

0558A 

0570A 

H 

D 

000576 

0428S 

0535M 

054 1M 

0543M 

0581 

0584 

HO 

D 

/COEF/ 

000033 

0428S 

0434S 

0535 

0630 

I 

I 

005001 

0438M 

0439 

0440 

0466M 

0467 

0490M 

0491 

0518M 

0519 

0551M 

0552 

0553 

0594M 

0595 

0596 

0617M 

0619 

INTERM 

I 

EXTERNAL 

000000 

0442 

I  SUM 

I 

005005 

0443M 

0448 

0484 

0526 

0529 

0539 

0548 

0559A 

0571A 

0579 

0589 

0599 

0623 

0634 

J 

I 

005006 

0465M 

0467 

0489M 

0491 

051 7M, 

0519 

0618M 

0619 

JL 

I 

005007 

0569M 

0570A 

0571 A 

0584 

0585 

JS 

I 

005010 

0545M 

0562M 

0569 

0580 

0581 

0582 

0584 

0585 

KEY 

I 

ARGUMENT 

000004 

0423S 

0588 

0616 

0620 

0645 

KR 

I 

ARGUMENT 

000003 

0423S 

0441A 

0485 

0486 

0487 

0488 

0491 

0493 

0494 

0495 

0496 

0498 

0499 

0500 

0501 

0507 

0508 

0537 

0632 

0637 

0647 

0648 

0649 

0650 

0651 

0652 

0653 

KS 

I 

005011 

0448M 

0485 

0486 

0487 

0488 

0491 

0493 

0494 

0495 

0496 

0498 

0499 

0500 

0501 

0507 

0508 

0537 

0632 

0637 

L 

D 

/INTEL/ 

000024 

0428S 

0434S 

0474 

0480 

0510 

0528 

0604 

0615 

0646M 

0647 

0652 

L2Z 

D 

/INTEL/ 

000034 

0428S 

0434S 

0449 

0450 

0468 

0471 

0472 

0481 

0482 

0492 

0604 

M 

I 

005013 

0556M 

0557 

0563 

0570 

0571 

0572 

0573 

0574 

0575 

0585 

0587 

0591 

0592 

MM 

D 

/COEF/ 

000020 

0428S 

0434S 

0538 

0633 

MS 

D 

/COEF/ 

000024 

0428S 

0434S 

0539 

0634 

MU 

D 

/COEF/ 

000004 

0428S 

0434S 

0480 

0528 

0532 

0540 

0604 

0615 

0625 

0635 

0646 

0652 

MUP 

I 

005014 

0555M 

0556 

N 

I 

005015 

0533M 

0535 

0542M 

0543 

0577M 

0578 

0579 

0581 

0582 

0584 

0585 

0587 

0589 

0628M 

0630 

0641M 

0643 

NO 

I 

005016 

0572M 

0574M 

0576 

0577 

0597 

N5 

I 

005017 

0578M 

0579M 

0585 

0586 

NMU 

I 

005020 

0554M 

0555 

NN 

I 

005021 

0560M 

NTER 

I 

005022 

0526M 

0533 

0542 

0554 

0558A 

0559A 

0570A 

0571 A 

0576 

0577 

0627A 

0628 

0640A 

0641 

NTERM 

I 

/COEF/ 

000030 

0433S 

0434S 

0526 

Q 

D 

000716 

0428S 

0558A 

0570A 

0581 

0584 

0627A 

0630 

0640A 

0643 

Q1 

D 

005023 

0428S 

0544M 

0558A 

0570A 

0626M 

0627A 

0639M 

0640A 
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(0634)  IF  (ISUM.EQ.3)  Z7=MS 

(0635)  Z7=Z7*MU»Z8*Z8/R 

(0636)  Z9=C2*SINB/( 1+RH0) 

(0637)  Z9=SL(KR ,KS+ 1 )* (COSU-Z9*TAU )+SL(KR ,KS+2)* (SINU-Z9*SIGMA )+ 

(0638)  +  SL(KR,KS4-3)*SINB 

(0639)  Q1 =1 

(0640)  CALL  QUE(0,Q,Z9,NTER,Q1 , 1 ) 

(0641 )  DO  570  N=2 ,NTER 

(0642)  Z7=Z7*Z8 

(0643)  570  DELF2=DELF2-Z7*Q( 1  ,N) 

(0644)  600  CONTINUE 

(0645)  IF  (KEY .EQ . 1 )  GOTO  650 

(0646)  L=L+L**3/ (MU*MU )*DELF2 

(0647)  VAR(6,KR)=L 

(0648)  650  VEL( 1 ,KR)=-DE(2) 

(0649)  VEL(2,KR)=DE( 1 ) 

(0650)  VEL ( 3 ,KR ) =-DE ( 4 ) 

(0651)  VEL(4,KR)=DE(3) 

(0652)  VEL(5,KR)=MU«MU/L»»3-DE(5) 

(0653)  VEL(6,KR)=0 

(0654)  RETURN 

(0655)  END 


A 

D 

004715 

0428S 

0622 

0528M 

0531 

0532 

0537 

0540 

0586 

ALPHA 

D 

004721 

0428S 

0481M 

0486M 

0511 

0512 

0513 

0520M 

0567 

0568 

0607 

0610 

AYBEE 

R  EXTERNAL 

000000 

0559 

0571 

B 

D 

000006 

0428S 

0559A 

0571 A 

0582 

0584 

0585 

B 1 

D 

004725 

0428S 

0547M 

0548M 

0559A 

0571A 

BETA 

D 

004731 

0428S 

0482M 

0487M 

0514 

0515 

0516 

0521M 

0567 

0568 

0607 

0611 

C 

D 

000246 

0428S 

04  39M 

0453M 

0454M 

0455M 

0456M 

0457M 

0460M 

0461M 

0462M 

0463M 

0464M 

0467M 

0469M 

0470M 

0471M 

0472M 

0473M 

0475M 

0476M 

0477M 

0478M 

0479M 

0480M 

0491M 

0493M 

0494M 

0495M 

0496M 

0497M 

0498M 

0499M 

0500M 

0501M 

05^2M 

0505M 

0506M 

0507M 

0508M 

0509M 

051  1M 

0512M 

051  3M 

0514M 

051 5M 

0516M 

0519M 

0619 

Cl 

D  /INTEL/ 

000044 

0428S 

04  3HS 

0519 

0520 

0521 

C2 

D  /INTEL/ 

000050 

0428S 

0434S 

0467 

0481 

0482 

0488 

0505 

0506 

0507 

0508 

0509 

0621 

0636 

CM 

C 

004735 

0428S 

0549M 

0565 

0567M 

0589 

0595 

0597 

CM1 

D 

004741 

0428S 

0565M 

0567 

0568 

0591 

COSU 

D  /INTEL/ 

000064 

0428S 

0434S 

0621 

0637 

D 

D 

004745 

0428S 

0560M 

0564M 

0575M 

0581 

0584 

0587M 

D 1 

D 

004765 

0428S 

0546M 

0560 

0563M 

0564 

DE 

D 

000466 

0428S 

0652 

0  4  4  OM 

0619M 

0648 

0649 

0650 

0651 

DELF2 

D 

004771 

0428S 

0530M 

0589M 

0616M 

0630M 

0643M 

0646 

DF 

D 

000516 

0428S 

0553M 

0591M 

0592M 

0595M 

0610M 

061 1M 

0612M 

0614M 

0615M 

0619 
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(0581) 

(0582) 

(0583) 

(0584) 

(0585) 

(0586) 

(0587) 

(0588) 

(0589) 

(0590) 

(059D 

(0592) 

(0593) 

(0594) 

(0595) 

(0596) 

(0597) 

(0598) 

(0599) 

(0600) 

(0601) 

(0602) 

(0603) 

(0604) 

(0605) 

(0606) 

(0607) 

(0608) 

(0609) 

(0610) 

(0611) 

(0612) 

(0613) 

(0614) 

(0615) 

(0616) 

(0617) 

(0618) 

(0619) 

(0620) 

(0621) 

(0622) 

(0623) 

(0624) 

(0625) 

(0626) 

(0627) 

(0628) 

(0629) 

(0630) 

(0631) 

(0632) 


Z9=D*H(N)*Q(JS,N) 

Z8=Z9*B(JS,N) 

F0( 1 )=F0( 1 )+Z8 

FO(3)=FO(3)+D«H(N)«Q(JL,N)»B(JS,N) 

F0(4)=F0(4)+(M-N5- 1 )/X*Z9* (B( JS,N)+B( JL,N) ) 

F0(5)=F0(5)+N5/A«Z8 

D=-D*(N-M+1 )/(N+M+2) 

IF  (KEY.EQ.1)  GOTO  400 

IF  (N.GT.2.0R .ISUM.GT . 1 )  DELF2=DELF2+Z8*CM 
400  CONTINUE 

DF(1)=DF(1)+F0(1 )*M*CM1 
DF(2)=DF(2)-F0( 1 )*M*SM1 
FO ( 1 ) = 0 
DO  410  1=3,5 
DF(I)=DF(I)+FO(I)*CM 
410  F0(I)=0 

IF(CM»CM+SM»SM.LT.1 .OD-8.AND.NO.GT.2)GOTO  510 
500  CONTINUE 

510  IF  (ISUM.NE. 1 )GOTO  530 
C 

C  SECOND  ORDER  CONTRIBUTION  OF  SECOND  HARMONIC 
C 

RH02=RH0»RH0 

Z9= (MU*GAM2P)**2*L2Z/ (48*L**3) 

Z8=5*ZETA«(URH02«(-3.6+RH02)) 

Z7=4»( URH02«(-6+9»RH02) ) 

Z6=ALPHA*ALPHA-BETA*BETA 

Z5=2-30»RHO2 

V=Z9* ( (Z8+Z7)*ZETA+5* (- 1+RH02* (2+7*RH02) )-Z6#Z5) 

DF( 1 )=DF( 1 )-2*Z9*Z5* ALPHA 
DF(2)=DF(2)+2«79*Z5*BETA 

DF(3)=DF(3)+Z9*4»PH0«(5*(1+7»RH02)+ZETA»(12«(-U3*RH02)+ZETA» 
+  (-9+5*RH02) )+15*Z6) 

DF ( 4)=DF ( 4)-ZETA* (ZETA*Z9* ( 2*Z8+Z7)-7*V ) 

DF(5)=DF(5)-5*MU*V/ (L*L) 

IF  (KEY.NE.1)  DELF2=DELF2+V 
530  DO  540  1=1,5 
DO  540  J=1 ,5 

540  DE ( I ) =DE ( I ) +DF (J)*C(J,I) 

IF  (KEY.EQ.1)  GOTO  600 
SINB=C2* (SIGMA*SINU+TAU*COSU ) 

R=A»ZETA»ZETA/ (ECOSF+ 1 ) 

IF  (ISUM.NE. 1)  GOTO  560 
Z8=RE/R 

Z7=-MU/R*Z8*Z8 
Q 1  =  1 

CALL  QUE(0,Q,SINB,NTER,Q1 , 1 ) 

DO  550  N=3,NTER 
Z7=Z7*Z8 

550  DELF2=DELF2-Z7*H0(N)*Q( 1 ,N) 

GOTO  600 

560  Z8=R/SL(KR ,KS+4 ) 
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(0528) 

A=L*L/MU 

(0529) 

IF  (ISUM.NE. 1 )G0T0  130 

(0530) 

DELF2=0 

(053D 

Z9=R£/A 

(0532) 

Z8=MU/A*Z9 

(0533) 

DO  120  N=2,NTER 

(0534) 

Z8=Z8*Z9 

(0535) 

120 

H(N)=-H0(N)*Z8 

(0536) 

GOTO  150 

(0537) 

130 

Z9=A/SL (KR ,KS+  4 ) 

(0538) 

Z7=MM 

(0539) 

IF ( ISUM . EQ . 3 ) Z7=MS 

(0540) 

Z8=MU»Z7/A»Z9 

(0541) 

H(1)=Z8«Z9 

(0542) 

DO  140  N=2,NTER 

(0543) 

140 

H(N)=H(N-1)»Z9 

(0544) 

150 

Q 1  =  1 

(0545) 

JS=1 

(0546) 

D1  =2 

(0547) 

B 1  =  1 

(0548) 

IF(ISUM.EQ.1)B1sX 

(0549) 

CM  =  1 

(0550) 

SM=0 

(0551) 

DO  160  1=1,6 

(0552) 

F0(I)=0 

(0553) 

160 

DF ( I ) =0 

(0554) 

NMU=NTER+1 

(0555) 

DO  500  MUP= 1 ,NMU 

(0556) 

M=MUP- 1 

(0557) 

IF  (M.GT.O)GOTO  180 

(0558) 

CALL  QU£(0, Q, GAMMA, NTER, Q1, 1 ) 

(0559) 

CALL  AYBEE(0,B,X,NTEF,B1 , 1 ,ISUM) 

(0560) 

D=D1/2 

(0561 ) 

GOTO  200 

(0562) 

180 

JS  =  3-vJS 

(0563) 

D1=D1/(2»M) 

(0564) 

D  =  D1 

(0565) 

CM1 =CM 

(0566) 

SM1=SM 

(0567) 

CM=CM1»ALPHA-SM1»BETA 

(0568) 

SM=SM1*ALPHA+CM1*BETA 

(0569) 

rvj 

o 

o 

JL=3-JS 

(0570) 

CALL  QUE(M*1 ,Q, GAMMA, NTER, Q1 ,JL) 

(0571  ) 

CALL  AYBEE (M+ 1,B,X,NTER,B1,JL, ISUM ) 

(0572) 

N0=M 

(0573) 

IF  (M.GE.2)  GOTO  300 

(0574) 

N0=M+2 

(0575) 

D=-D/ (2*M+2) 

(0576  ) 

300 

IF (N0.GT.NTER )GOTO  500 

(0577) 

DO  400  N=N0,NTER,2 

(0578) 

N5=N 

(0579) 

IF ( ISUM.EQ. 1 )N5=-N-1 
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(0475) 

(0476) 

(0477) 

(0478) 

(0479) 

(0480) 

(0481) 

(0482) 

(0483) 

(0484) 

(0485) 

(0486) 

(0487) 

(0488) 

(0489) 

(0490) 

(0491) 

(0492) 

(0493) 

(0494) 

(0495) 

(0496) 

(0497) 

(0498) 

(0499) 

(0500) 

(0501) 

(0502) 

(0503) 

(0504) 

(0505) 

(0506) 

(0507) 

(0508) 

(0509) 

(0510) 

(051D 

(0512) 

(0513) 

(0514) 

(0515) 

(0516) 

(0517) 

(0518) 

(0519) 

(0520) 

(0521) 

(0522) 

(0523) 

(0524) 

(0525) 

(0526) 

(0527) 


C(4,2)=-XI*Z9 
C(4, 1 )=-ETA*Z9 
C(4,4)=0 
C(4,3)=0 

C(4,5)=(1-ZETA)»Z9 
C(5,5)=2*L/MU 
ALPHA =C2«V1»L2Z 
BETA=C2*V2*L2Z 
GAMMA =RHO 

IF ( ISUM . EQ . 1 )G0T0  100 

Z7=SL(KR ,KS+ 1 )*TAU+SL(KR ,KS+2 )*SIGMA 

ALPHA =SL(KR  ,KS+3 )*ALPHA+SL(KR  ,KS+ 1  )*XI-SL(KR  ,KS+2)*ETA-V1 *Z7 
BETA=SL(KR ,KS+3 )*BETA-SL(KR ,KS+ 1 )*ETA-SL(KR ,KS+2)"XI-V2»Z7 
GAMMA =SL  OCR  ,KS+3 )*GAMMA-C2*Z7 
00  90  J=  1 ,5 
DO  90  1=1,3 

90  C(I,J)=SL(KR,KS+3)*C(I,J) 

Z7=Z7/L2Z 

C(1,2)=C(1 ,2)+SL(KR ,KS+ 1 )-Z7* (TAU+2*V 1*XI ) 

C(  1 , 1 )=C( 1 , 1 )-SL(KR ,KS+2)+Z7* (SIGMA-2* VI *ETA) 

C ( 1 , 4)=C ( 1 , 4)-SL(KR ,KS+2)*V 1+ETA*Z7 
C ( 1 , 3) =C ( 1 , 3)-SL(KR ,KS+ 1 )*V 1-XI*Z7 
C(1,5)=C(1 ,5)+2*V1*Z7 

C (2,2)=C(2,2)-SL(KR ,KS+2)+Z7* (SIGMA-2*V2*XI ) 

C(2, 1 )=C(2, 1 )-SL(KR ,KS+ 1  HZ7" (TAU-2*V2*ETA) 
C(2,4)=C(2,4)-SL(KR ,KS+2)*V2+XI*Z7 
C(2,3)=C(2, 3)-SL(KR ,KS+1 )*V2+ETA*Z7 
C(2,5)=C(2,5)+2*V2*Z7 
Z9=Z7/(URHO) 

Z8=2*RHO»Z9 

C(3,2)=C(3,2)-C2*XI*Z8 

C(3,1)=C(3,1 )-C2*ETA*Z8 

C ( 3, 4)=C ( 3 ,4)-C2* (SL(KR ,KS+2)-SIGMA*Z9) 

C ( 3, 3)=C ( 3, 3)-C2* (SL(KR ,KS+1 )-TAU*Z9) 
C(3,5)=C(3,5)+C2*2*RHO*Z9 
100  Z9=1/(2*L*(HZETA)) 

CO  ,2)=C(1 ,2)-ALPHA*XI*Z9 
C( 1 , 1 )=C( 1 , 1 )-ALPHA*ETA*Z9 
C(1,5)=C(1,5)-2»ALPHA«ZETA*Z9 
C(2,2)=C(2,2)-BETA*XI»Z9 
C(2, 1 )=C(2, 1 )-BETA*ETA*Z9 
C(2,5)=C(2,5)-2*BETA»ZETA*Z9 
DO  110  J=1,5 
DO  110  1=1,2 
110  C(I,J)=C1*C(I,J) 

ALPHA=C1*ALPHA 

BETA=C1*BETA 

C 

C  RECURSIVE  CALCULATION  OF  DERIVATIVES  OF  PERTURBING 
C  POTENTIAL  WITH  RESPECT  TO  XI , ETA , SIGMA  ,TAU  ,L 
C 

NTER  =NTERM ( ISUM ) 

X= 1/ZETA 
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(0422)  C 

(0423) 

(0424)  C 

(0425)  C 

(0426)  C 

(0427)  C 

(0428) 

(0429) 

(0430) 

(0431  ) 

(0432) 

(0433) 

(0434) 

(0435) 

(0436) 

(0437)  C 

(0438) 

(0439) 

(0440)  70 

(0441  ) 

(0442) 

(0443) 

(0444)  C 

(0445)  C 

(0446)  C 

(0447)  C 

(0448) 

(0449) 

(0450) 

(0451) 

(0452) 

(0453) 

(0454) 

(0455) 

(0456) 

(0457) 

(0458) 

(0459) 

(0460) 

(0461  ) 

(0462) 

(0463) 

(0464) 

(0465) 

(0466) 

(0467)  80 

(0468) 

(0469) 

(0470) 

(0471 ) 

(0472) 

(0473) 

(0474) 


SUBROUTINE  DERIV (KR .KEY ) 

KEY=1:  DO  NOT  CALCULATE  SECOND  ORDER  CORRECTIONS  TO  L 
KEY=2:  DO 

DOUBLE  PRECISION  VAR ( 6 , 1 7) , VEL ( 6 , 1 7 ) ,XI ,ETA .SIGMA ,TAU ,L .LAMBDA , 

+  ZETA,L2Z,RH0,C1,C2,V1 , V2, Z5 , Z6 ,Z7 , Z8.Z9, ALPHA, BETA. GAMMA, C (6,6), 
+  DRAD ,MU .GAM2N , RE ,MM ,MS ,X ,A , B1 ,D,D1,Q1 ,CM1 ,SM1 ,CM,SM,R,B(2,20) , 

«■  Q(2,2O),F0(6),DF(6),DE(6),SL(17,8),H(2O),H0(2O),TT,DELF2,SINB, 

♦  ECLON , U ,SINU ,COSU , ECOSF , ESINF , EC ,ES .GAM2P ,RH02 , V 
DIMENSION  NTERM( 3) 

COMMON/INTEL/ETA, XI, TAU, SIGMA, LAMBDA, L.ZETA.L2Z.RH0, Cl ,C2, 

♦  ECLON  ,U  ,COSU , SINU  .ECOSF  .ESINF , EC  ,ES  .GAM2P 

+  // VAR , VEL/COEF/DRAD ,MU .GAM2N , RE ,MM ,MS .NTERM .HO/LUNSOL/TT ,SL 

DO  70  1=1,5 
C (5 , I) =0 
DE ( I ) =0 
CALL  XFR(KR ,2) 

CALL  INTERM 

DO  600  ISUM= 1 , 3  /*  UZONALS,  2=MOON,  3=SUN 

JACOBIAN  OF  ALPHA, BETA, GAMMA, X, A  WITH  RESPECT  TO 
ETA, XI, TAU, SIGMA ,L 

KS=4*(ISUM-2) 

V1=(XI*TAU-ETA»SIGMA)/L2Z 

V2=(-XI«SIGMA-ETA»TAU)/L2Z 

Z8=V1/(URHO) 

Z9=2«RHO«Z8 

C ( 1 ,2) =TAU+Z9*XI 

C(  1 ,  1  )=-SIGMA>-Z9«ETA 

C( 1 ,4)=-ETA-SIGMA»Z8 

C(1,3)=XI-TAU*Z8 

C( 1 ,5)=-Z9 

Z8  =  V2/(  URHO) 

Z9^2*RHO«Z8 
C(2,2)=-SIGMA+Z9*XI 
C (2, 1 )=-TAU+Z9*ETA 
C ( 2, 4) =-XI-SIGMA*Z8 
C(2,3)=-ETA-TAU*Z8 
C (2,5) =-Z9 
DO  80  J=1 ,5 
DO  80  1=1,2 
C  ( I ,  J  )  =C2*C ( I , J ) 

Z9=2* ( 1-RHO )/L2Z 
C(3,2)=-XI*Z9 
C ( 3 , 1 )=-ETA*Z9 
C ( 3 , 4 )=-2*S IGMA/L2Z 
C ( 3, 3) =-2*TAU/L2Z 
C(3,5)=Z9 

Z9=-1/(L»ZETA»ZETA) 
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D  ARGUMENT 

000007 

0382S 

0384S 

0390 

0391 

0392M 

0410 

0412 

0419M 

D  ARGUMENT 

000004 

0382S 

0384S 

0394M 

0395M 

0397M 

0399 

0407 

0408M 

0410M 

0412M 

0415M 

I  ARGUMENT 

00001 1 

0382S 

0386 

I  ARGUMENT 

000010 

0382S 

0394 

0395 

0397 

0399 

0407 

0408 

0410 

0412 

0415 

I  ARGUMENT 

000003 

0382S 

0391 

0393 

0394 

0395 

0397 

0406 

0407 

0408 

0409 

0410 

0411 

0412 

0413 

0415 

I 

000371 

0393M 

0397 

0399 

0400M 

0401 

041 3M 

0415 

0417M 

0418 

I  ARGUMENT 

000006 

0382S 

0401 

0409 

0411 

0418 

D 

000372 

0384S 

0390M 

0394 

0397 

0398M 

D 

000376 

0384S 

0391M 

0392 

0395 

0397 

0398 

0399 

D  ARGUMENT 

000005 

0382S 

0384S 

0385 

D 

000402 

0384S 

0385M 

0397 

0412 

0415 

0419 

000106 

0397D 

0402 

000172 

0396 

0401D 

000177 

0386 

0406D 

000227 

0406 

0410D 

000263 

0415D 

0418 

000344 

0414 

0418D 
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(0422)  C 
(0423) 
(0424)  C 
(0425)  C 
(0426)  C 
(0427)  C 
(0428) 
(0429) 
(0430) 

(0431  ) 

(0432) 

(0433) 

(0434) 

(0435) 

(0436) 

(0437)  C 

(0438) 

(0439) 

(0440)  70 

(0441  ) 

(0442) 

(0443) 

(0444)  C 

(0445)  C 

(0446)  C 

(0447)  C 

(0448) 

(0449) 

(0450) 

(0451) 

(0452) 

(0453) 

(0454) 

(0455) 

(0456) 

(0457) 

(0458) 

(0459) 

(0460) 

(0461  ) 

(0462) 

(0463) 

(0464) 

(0465) 

(0466) 

(0467)  80 

(0468) 

(0469) 

(0470) 

(0471) 

(0472) 

(0473) 

(0474) 


SUBROUTINE  DERIV (KR ,KEY  ) 

KEY= 1 :  DO  NOT  CALCULATE  SECOND  ORDER  CORRECTIONS  TO  L 
KEY=2:  DO 

DOUBLE  PRECISION  VAR ( 6 , 1 7 ) , VEL ( 6 , 1 7 ) ,XI ,ETA , SIGMA ,TAU ,L .LAMBDA , 

+  ZETA,L2Z,RH0,C1 ,C2,V1 , V2, Z5 ,Z6 ,Z7 ,Z8, Z9 , ALPHA ,BETA ,GAMMA ,C (6 ,6) , 
+  DRAD ,MU ,GAM2N ,RE ,MM ,MS ,X , A ,B1 ,D,D1 ,Q1 ,CM1 ,SM1 ,CM,SM,R,B(2,20) , 

+  Q(2,20),F0(6),DF(6),DE(6),SL(17,8),H(20),H0(20),TT,DELF2,SINB, 

♦  ECLON ,U ,SINU ,COSU ,ECOSF ,ESINF , EC ,ES ,GAM2P , RH02 , V 
DIMENSION  NTERM(3) 

COMMON/INTEL/ETA, XI, TAU, SIGMA, LAMBDA, L,ZETA,L2Z,RH0, Cl ,C2, 

+  ECLON ,U ,COSU ,SINU ,ECOSF ,ESINF,EC ,ES ,GAM2P 

+  //VAR , VEL/COEF/DRAD ,MU ,GAM2N ,RE ,MM ,MS ,NTERM ,HO/LUNSOL/TT ,SL 

DO  70  1=1,5 

C(5,I)=0 

DE(I)=0 

CALL  XFR (KR , 2 ) 

CALL  INTERM 

DO  600  ISUM= 1 , 3  /»  1 =ZONALS ,  2=M00N,  3=SUN 

JACOBIAN  OF  ALPHA, BETA, GAMMA, X, A  WITH  RESPECT  TO 
ETA, XI, TAU, SIGMA, L 

KS=4* (ISUM-2) 

V 1  =  (XI*TAU-ETA*SIGMA ) /L2Z 
V2= (-XI*SIGMA-ETA*TAU )/L2Z 
Z8=V1/( 1+RHO) 

Z9=2«RHO*Z8 
C ( 1 , 2) =TAU  +  Z9*XI 
C  ( 1 , 1 )  =  -SIGMA  +  Z9*ETA 
C(  1 ,4)=-ETA-SIGMA»Z8 
C(  1 ,3)=XI-TAU*Z8 
C ( 1 , 5) =-Z9 
Z8  =  V2/(URHO) 

Z9=2*RH0*Z8 
C(2,2)=-SIGMA+Z9*XI 
C (2, 1 )=-TAU+Z9*ETA 
C ( 2, 4) =-XI-SIGMA*Z8 
C ( 2 , 3 ) = -ET  A -T  AU  *  Z8 
C(2,5)=-Z9 
DO  80  J=1 ,5 
DO  80  1=1,2 
C ( I , J) =C2*C ( I , J ) 

Z9=2»(1-RHO)/L2Z 
C ( 3*2 ) =-XI*Z9 
C ( 3  *  1 ) =-ETA*Z9 
C ( 3 , 4) =-2*SIGMA/L2Z 
C ( 3  *  3 ) =-2*TAU/L2Z 
C(3,5)=Z9 

Z9  =  —  1 / (L"ZETA*ZETA ) 


A-22 
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QUE 


Calculates  the  quantities  Q  required  for  the  harmonics. 


AYBEE 


Calculates  the  quantities  A  and  B  required  for  the  harmonics. 


DERIV  (KEY) 


Loads  the  mean  elements  with  XFR  and  INTERM;  and  in  separate  calculations  for  zonals, 
moon,  and  sun: 

1.  Finds  the  Jacobian  matrix  of  the  variables  ALPHA,  BETA,  GAMMA,  X,  A  (in  terms  of 
which  the  derivatives  of  the  Hamiltonian  are  naturally  determined)  with  respect  to  the 
mean  elements. 

2.  Calculates  the  derivatives  of  the  perturbing  potential  as  well  as  the  potential  value  itself. 
The  second  order  of  the  second  harmonic  is  included. 

3.  If  KEY=2  (at  entry  to  the  program  only),  finds  the  instantaneous  value  of  the  perturbing 
potential  and  uses  it  to  calculate  the  remaining  part  of  the  second-order-short-period  cor¬ 
rection  for  L.  If  KEY=1.  this  step  is  omitted. 

4.  Use  the  Jacobian  and  the  derivatives  to  accumulate  the  derivatives  of  the  mean  elements 
required  for  the  next  integration  step,  storing  these  in  VEL. 


EXTRAP  (KROW1,  DROW) 

Extends  the  difference  table  DEL  from  the  preceding  row  KROW1  to  the  current  row 
KROW,  calculates  the  step  in  VAR.  and  calls  DERIV  for  the  new  row. 


NSWC  TR  84-25 


SUBROUTINE  TAXONOMY 


NAME 

CALLS 

IS  CALLED  BY 

GPS 

ENTER 

(MAIN) 

OSTCOR 

ELTRAN 

READ 

DERIV 

EXTRAP 

ENTER 

INTERM 

XFR 

CiPS 

OSTCOR 

INTERM 

XFR 

GPS 

ELTRAN 

INTERM 

XFR 

GPS 

INTERM 

ENTER 

OSTCOR 

ELTRAN 

XFR 

ENTER 

OSTCOR 

ELTRAN 

READ 

(iPS 

OUT 

DERIV 

AYBI  E 

DERIV 

DERIV 

INTERM 

GPS 

XFR 

OUE 

AY  BEE 

EXTRAP 

EXT RAP 

DERIV 

(.PS 

USES  COMMON 

(BLANK) 

COEF 

COREL 


COEF 

COREL 

INTEL 

COEF 

COREL 

(BLANK) 

COEF 

COREL 

INTEL 

COEF 

INTEL 


(BLANK) 

COREL 

INTEL 

LUNSOL 


COEF 

INTEL 

LUNSOL 


(BLANK) 

COEF 
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VARIABLES 

AND  COMMON  BLOCKS 

(BLANK) 

Integration  quantities 

VAR  (6,17) 

The  six  integration  variables 
(mean  elements)  for  17  con¬ 
secutive  time  steps 

VEL  (6,17) 

The  derivatives  in  a  similar 
array 

DEL  (18,6,17) 

The  differences  used  in  the 
numerical  integration 

NEQ 

Number  of  equations 

KMAX 

Order  of  numerical  integration 
procedure 

COEF 

Numerical  coefficients 

DRAD 

Degrees/radians 

MU 

Gravitational  const 

GAM2N 

1 2*J2*(MU*RE)**2 

RE 

Earth  radius 

MM 

Lunar  mass 

MS 

Solar  mass 

NTERM  (3) 

Highest  order  harmonic: 
zonal,  lunar,  solar 

HAR(20) 

Coefficients  of  zonal  harmonics 

ADAM  (17) 

Adams  coefficients 

COREL 

Coordinates  and  elements 

X,Y,Z,VX,VY,VZ  (or  C(6)) 

Coordinates 

A,ECC,INC,ELL,G,H  (or  E(6)) 

Elliptic  elements 

0(6) 

Osculating  elements 

INTEL 

Intermediate  quantities  in  element  transformations 

ETA, XI, TAU, SIGMA, LAMDA,L  (or  W(6)) 

Working  elements  (osculating 
on  mean) 
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114  variables,  such  as  ZETA.  RHO.  U.  Intermediate  variables 

ECOSF-  see  listing) 

LUNSOL  Lunar  and  solar  coordinates 

TSL  Time 

SL  (17,  8)  For  17  time  steps,  the  three 

components  of  the  unit  vector 
to  the  moon  and  the  lunar 
distance;  followed  by  the 
corresponding  4  quantities 
for  the  sum 

The  variables  not  found  in  common  declarations  are  few  in  number  and  usually  nonce  variables. 

An  exception  is  the  variable  KS  that  counts  the  .steps  in  the  integration.  The  variable  TIM 

(which  is  defined  but  not  used)  is  the  time;  its  definition  is 

TIM=STEP*(KS-2)+TIMO  (seconds) 

External  procedures  not  in  the  subroutine  list  (system  subroutines  or  functions): 


EXIT 

DSQRT 

MOD 

DABS 

SQRT 

DATAN2 

DCOS 

DMOD 

DSIN 

DSIGN 

All  of  the  DOUBLE  PRECISION  declarations  and  functions  are  to  produce  a  64-bit  value 
(48-bit  mantissa).  On  a  machine  with  sufficient  word  length,  all  the  DOUBLE  PRECISION 
declarations  could  be  replaced  with  REAL  and  all  the  double  functions  with  their  single¬ 
precision  counterparts. 
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TRANSFORMATION  FROM  CARTESIAN  COORDINATES  TO  OSCULATING  CANONICAL 
VARIABLES 


[This  is  merely  summarized  here,  since  most  of  it  is  to  be  found  in  any  text  on  orbital 
mechanics  and  the  rest  is  obvious  from  the  definition  of  the  canonical  variables.  In  the  program, 
it  is  found  in  subroutine  ENTER  (lines  147-178)].  Vectors  are  defined  by  an  underline 
unit  vectors  by  a  caret  A  . 


Angular  momentum  h  =  r  x  v 

magnitude  h  =  \/h^  +  h*  +  h* 


Unit  vector 


Hence, 


Unit  vector 


Unit  vector 


A 

h  = 


sin  i  sin  h 
-sin  i  cos  h 
cos  i 


A 

L, 


-  i|/u+p)  \ 

-  Is<c«l+0)J 

/ 


l2 


-I-1C/(1+P)\ 
-I*/(I+P)  j 


A  A  A 

L,.  L2.  and  h  are  mutually  orthogonal  and  define  a  Cartesian  system  rotated  from  the  XYZ 

A  A 

system  by  the  angle  i  about  the  line  of  nodes.  Therefore,  L,  and  L2  define  the  orbital  plane, 

A 

the  radius  vector  makes  the  angle  u  with  L,,  so  that 


A  A 

cos  u  =  r«L, 


A  A 

sin  u  =  r-L2 


PREVIOUS  PAGE 
IS  BLANK 
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[The  formulation  here  differs  from  the  customary  one  in  terms  of  the  basis  defined  by  N  (the 

A  a  A  A  A  A 

vector  along  the  line  ot  nodes),  h,  and  h  x  N.  When  the  inclination  is  small,  N  and  h  x  N  are 

A  A 

ill-detined,  a  disadvantage  not  suffered  by  1^  and  L2.  The  change  entails  also  using  the  long¬ 
itudes  (u,  a,  X)  instead  of  the  anomalies  (f,  H,  12)]. 

From  here  on  the  algorithm  is  straightforward: 

h  =  yv  -  zv  \ 

*  *  y  I 


h  =  zv  -  xv 

y  x  z 


h  =  xv  -  yv 

z  y  x 


h  =  \/  h2  +  h2  +  h2 

x  y  z 


L  =  -  h  /h 

S  x ' 


L  =  -  h  /h 

C  y ' 


P  =  h  /h 


x2  +  y2  +  z2 


cos  u  =  (x  +  zls/(l+p))/r 


sin  u  =  (y  +  zlt. /( l+p))/r 


e  cos  f  =  (h2/pr)  -  1 


e  sin  f  =  (h/pr)(xvx  +  yv  +  zvz ) 


■r.  •  .  '  *  V  .  . 


•  ■  _  »  _  •  M  * 
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e2  =  (e  sin  f)2  +  (e  cos  f)2 
f  = 

e  cos  cj  =  e  cos  f  cos  u  +  e  sin  f  sin  u 
e  sin  cj  =  e  cos  f  sin  u  -  e  sin  f  cos  u 

h2  ^  ^ 

—  cos  a  =  cos  u  +  e  cos  cj  +  e  sin  co.e  sin  f/(l  +  f) 

Mr 

h2  ^ 

—  sin  a  =  sin  u  +  e  sin  cj  -  e  cos  cj  .e  sin  f/(l  +  f) 

Mr 

,  (h2/pr)  sin  a 

a  =  tan" 1  — — - 

(Ir/pr)  cos  a 

L  =  h  /? 

X  =  a  -  e  cos  cj  sin  a  +  e  sin  cj  cos  a 

£  =  V  2L/(  1+f)  e  cos  cj 

r ?  =  -  V2L/(l+f)  e  sin  cj 

a  =  V2Lf/(l+p)  Ic 

t  =  V2Lf/(I+p)  Ig 

TRANSFORMATION  FROM  ELLIPTIC  ELEMENTS  TO  OSCULATING  CANONICAL 
VARIABLES 

[This  is  the  first  section  of  the  ENTER  subroutine  (lines  134-145)  invoked  when  the  KEY 
argument  is  2.) 

L  =  -s/pa 


X  =  £  +  g  +  h 
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£  =  \/2L(l  -  f)  cos  (g  +  h) 

r?  =  -  V 2L  (1  -  f)  sin  (g  +  h) 
a  =  2  \/Lf  sin  i/2  cos  h 

r  =  -  2  v/Tf- sin  i/2  sin  h 


TRANSFORMATION  FROM  OSCULATING  CANONICAL  VARIABLES  TO 
ELLIPTIC  ELEMENTS 

[This  is  in  subroutine  OSTCOR  (lines  194-203).]  See  the  last  section  in  this  Appendix  for 
the  calculation  of  the  intermediate  quantities  used  below. 


a  =  L2/p 

e  =  Vi  -  r2 

i  =  cos’ 1  p 
h  =  tan-1  (-  r/a) 
g  +  h  =  tan'1  (-T ?/£) 

6  =  X  -g  -h 

TRANSFORMATION  FROM  OSCULATING  CANONICAL  VARIABLES  TO 
CARTESIAN  COORDINATES 

| This  is  subroutine  OSTCOR  (lines  205-212).]  Again  the  intermediate  quantities  described  in 
the  next  section  are  used. 


r  =  (1  -  e  cos  oj  cos  a  -  e  sin  cc  sin  a)L2/p 
x  =  r((l  -  r2/2Lf)  cos  u  -(ar/2Lf)  sin  u) 
y  =  r(~(o  t/2 Lf)  cos  u  +  (1  -  a2  /2Lf)  sin  u) 


i 

C-b 
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z  =  r  (a  sin  u  +  t  cos  u)  V ( l+p)/2Lf 
vx  =  ((a  -  X)  x  +  f  (-  oz  V'(l+p)/2Lf  -  py))L/r2 
v  =  ((a  -  X)  y  +  £  (px  +  rz  x/0+p)2Lf))L/r2 
vz  =  ((a  -  X)  z  +  f  (-  ry  +  ox)  V  (1  +  p)2Lf)L/r2 

CALCULATION  OF  INTERMEDIATE  QUANTITIES  FROM  CANONICAL  VARIABLES 
(This  is  the  subroutine  INTERM.) 

f  =  i  -  a2  +  v2  >/2L 

p  =  1  -  (a2  +  r2) 2Lf 
e  cos  cj  =  £  V(1  +  ?)/2L 

esincj  =  -r?V(l+  £)/2L 

Solve  for  a,  by  Newton-Raphson,  the  Kepler  equation: 

X  =  a  -  e  cos  a>  sin  a  +  e  sin  oj  cos  a 

u  =  a  +  2  tan-1  ((a  -  X)/(l  +  £  -  e  cos  tu  cos  a  -  e  sin  co  sin  a)) 
e  sin  f  =  e  cos  cu  sin  u  -  e  sin  c5  cos  u 
e  cos  f  =  e  cos  gj  cos  u  +  e  sin  to  sin  u 

c  =  vrrrjj/jL 

C2  =  Vo  +  P)2L£ 

V,  =  J2  <n  Re)2/(2L?)4 


C-7 


-A. -»  J  ^  -*■  - *  •*.  -U 


.  -  %“  V*  %"  *, 


J  -*  ^  A 
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APPENDIX  D 

DIFFERENTIAL  RELATIONSHIPS  AND  THE  TRANSFORMATION 
BETWEEN  OSCULATING  AND  MEAN  CANONICAL  VARIABLES 
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DIFFERENTIAL  RELATIONSHIPS  AND  THE  TRANSFORMATION 
BETWEEN  OSCULATING  AND  MEAN  CANONICAL  VARIABLES 


The  first-order  determining  function,  whose  partial  derivatives  give  the  transformation  between 
mean  and  osculating  variables,  is  given  by 


S,  =  6  (pRe)2  J2/(2Lf)3  j(-  1+3  P2)U,  + 

+  (  (1  +  P)/4Lf)(a2  -  r2)  Sup  +  2  or  Cup)[ 


Uv  =  f  -  £  +  e  sin  f 

S  =  sin  2u  +  e  sin  (u  +  cj)  +  (e/3)  sin  (3u  -  cj) 


Cup  =  cos  2u  +  e  cos  (u  +  cj)  +  (e/3)  cos  (3u  -  cj) 

The  advantage  of  this  dissection  is  that  none  of  the  last  three  quantities  defined  depend  upon 
a  or  r.  Quantities  such  as  e,  f.  and  C  are  to  be  considered  as  expressed  in  terms  of  the  canon¬ 
ical  variables  being  used.  Hence,  it  is  necessary  to  derive  a  miscellany  of  differential  relationships 
that  will  be  needed  in  the  process. 


Since  f2 


de  =  -  fd  f/e  = 


fe 


i+r 


2(£  d£  +  t?  dp)  dL 

s2  +  p2  T 


Also,  dw  =  d(tan_l 


(-  Hit))  = 


t?dg  -  l-dp 

$2  +  n2 


From  the  usual  Keplerian  orbit  relationships: 


8  =  E  -  e  sin  E 
r/'a  sin  f  =  f  sin  E 
r/a  cos  f  =  cos  E  -  e 

f 


D-3 
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Differentiation,  substitution,  and  rearrangement  yields 


if  =  (Z2/f3)  dC  +  ((Z  +  1)  sin  f)/f2)de 
vhere  Z  =  f2  a/r  =  1  +  e  eos  f. 
t  will  also  be  convenient  to  use 


ec  =  VTI  +  f)/2L  *  =  e,  * 

=  -  c,  rj 

rom  which  by  differentiation: 

dec  =  ct  d£  +  ec/(2L(l  +  f))  (-  2  ?  dL  -  £  d£  -  r?  dr?) 
des  =  -  c,  dr?  +  es/(2L(l  +  ?))  (~  2  f  dL  -  £  d£  -  7?  dr?) 
f  the  foregoing  are  applied, 

dUt  =  df(  1  +  e  cos  0  -  dC  +  sin  f  de 

=  <Z3/f3  -  1)  <dX  -  du)  +  (Z2  +  Z  +  r2)  sin  f  de/f2 


f2  +  n 2 


[<Z3/r3  -D(-  r?  d£  +  £  dr?) 


(Z2  +  Z  +  f2) 

+  2  - L- Le 


?(1  +  f) 


sin  f(£  d£  +  r?  dr?)] 


.  ,  2(Z2  +  Z  +  f2) 

+  (Z3/?3  -  I  )dX - —  e  sin  1  dL 

Lr  (1  +  f) 


he  coefficients  of  d£  and  dr?  are  not  convenient  for  computing,  since  at  low  eccentricities 
ley  become  the  ratios  of  small  quantities.  It  is  expedient  to  transform  them  using  the  easily 
erived  relations: 


-  r?  cos  f  +  £  sin  f  =  V  2L/<  1  +  f )  e  sin  u 
r?  sin  f  +  £  cos  f  =  \/  2L/(  1  +  f )  e  cos  u 

D-4 
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arriving  finally  at 


1  +  Z  +  Z2 

dU„  =  - -  12L(1  +  f)  sin  u  d£  +  cos  u  dr?) 


2Lf3 


2f+l  1  +  t  +  ?2 

e  sin  f  (£d£  +  7?d77) J  +  — — — -  (-  tj  d£  +  £  dr?) 


i+r 


2U 


(£  d£  +  r?  d7j)/L?  -  z:  e  sin  f  dL  +  (Z3/f3  -  l)dX 


Lf  (1  +  ft 


Next  Cup  may  be  written  =  cos  2u  +  ep(cos  u  +  (1/3)  cos  3u) 


+  e  (-  sin  u  +  (1/3)  sin  3u) 


so  that  dCup  =  -  (2  sin  2u  +  ec  (sin  u  +  sin  3u)  +  eg  (cos  u  -  cos  3u)  (df  +  d^ci 


+  (cos  u  +  (1/3)  cos  3  u)  dec  +  (-sin  u  +  (1/3)  sin  3u)des 
Note  that  2  sin  2  u  +  ec  (sin  u  +  sin  3u)  +  es  (cos  u  -  cos  3u)  =  2Z  sin  2  u; 


substitute,  rearrange,  and  make  the  same  transformation  as  for  U8.  The  result  is 


,  1  +  ? 

dCnn  =  -  (Z  sin  Z.i/Lf3)  [ (Z  +  1) -  (sin  u  d£  +  cos  u  dr?) 


u  p 


1  +  2f 


1  +  S 


e  sin  f(£  d£  +  r?  d7?))  +  ( 1  +  f  +  f2 )  (-  7?  d£  +  £  dr?)] 


+  c  | cos  u  +  (1/3)  cos  3  u)d£  +  (sin  u  -  (1/3)  sin  3  u)dr?] 


^  u  p  COS  _  U 


dL 


- — -  <£  d£  +  t?  dr?)+  7  -  — ; — —  |  2Z  (Z  +  1)  sin  2u.e  sin  f 


2L(  1  +  £) 


l?  (i  +  r> 


-  f2  (Cup  -  cos  2  u ) ]  -  2Z3  sin  2  u  dX/f3 


1 

'J 
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—  =  c,  (x'r  +  y'a)  2va/(2Lf)  +  c,c2z'  (-  2v2  p)/(  1  +  p)  -  2  0J/(2L(1  +  £)) 


~  =  (2  £/(2Lf))(-  c2  p  (x'r  +  y'o)  /(I  +  p)  - 


z'  (1  -  P)) 


87 

3rj 


-  (2i?/(2Lf))(-  c2  p(x'r  +  y'a)/(  1  +  p)  -  z'  (1  -  p)) 


87 

3a 


-  c2a  (x'r  +  y'a)/(2Lf  (1  +  p))  -  c2y'  -  2z'a/(2Lf) 


1  -  ,  > 

—  -  c2r  (x  r 

3r 


+  y'or)/(2Lf  (1  +  p))  -  c2x'  -  2z'r/(2Lf) 


87 

8L 


2c2p  (x'r  +  y'a)/(2Lf  (1  +  P))  +  2z'  (J  -  p)/(2Lf) 


9£ 


=  -  */L 


8f 

—  =  ( 1 
8L 


n/L 


8a 

8L 


2L  in 


The  other  partial  derivatives  in  the  matrix  are  zero. 
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f  =  1  -  (£2  +  172  )/2L, 

a 

=  L2/p 

we  may  write 

<*x  =  e,  (£  -  v,  r) 

a 

y 

=  ‘-'i  ("  1?  ~  v,  a) 

«2  =  c2  V,  (2Lf) 

dx  =  S',  (-  T}  -  v2  T) 

0y 

=  ‘-'i  (-  %  ~  v2  0) 

dz  =  c2  v2  (2Lr) 

yx  -  -  c2  T  7y  =  -  C2  a  7X  =  P 

and  easily  but  tediously  differentiate  these  to  find 


9a 

3?' 


9a 

9r? 


9a 


e,  x'-  c j  (r+2vif)(x'r  +  y'a)/(2U)  +  c1c2z'(r  +  2v]  £  p/(  1  +  p))  -  a£/(2L(l  +  f)) 


c,y' +  e,  (a  -  2v)  7?)(x'r  +  y'a)/(2Lf)  +  c,  c2z'  (-  0+2V,  17  p/(  1  +p))-  ai?/(2L(l  +f)) 


~  =  c,  (x  r  +  y'0)r?/(2Lf)-  c,  y'v,  +c,c2z'  (-  r?  -  ov,/(  1  +  p)) 


9  a 
dr 


da 


--  e,  (x'r  +  y'a)  £/(2L£)-  clX'v,  +  c,c2z'(£-  rvjil  +  p)) 


—  -  c,  (  x'r  +  y'o  )2v,  /( 2Lf )  +  e,  c2  z' (~  2v,  p)/(  1  +p)-  2af/(2L(l  +  f)) 


3d 


L'|>'  +  C1  2v2£)(\'r  +  y'a)/(2Lf)  +  CjC2z'(-  a  +  2v2  £p/(  1  +  p))  -  0$/(2L<  1  +  £)) 


3d 

3r? 


x'  +  e,  (r-  2v2  17 H x'r  +  y'o)/( 2Lf)  +  e, e2z’  (-  r+  2v2  77  p/C  1  +  p))  -  dr?/(2L(  1  +f)) 


-c,  (x'r  +  y'o)  £/(2LJ)-  c,  y'v2  +c,c2z'(-  £-  v2  a>/(  1  +p)) 


3d 
3o 

3d 

—  -  c(  (\  r  +  y  o)r?/(2Lf)-cIx'v2  +c,e2z'(-r?-  v2  r/<  1  +p)) 


G-4 


■3 

‘.<1 

■j 


.  «  •  *  H  -  »  "  k*4  '»***-*  •  '  "  •  "  *  ■*  «  ■  •  ■  »  •  ■ 


niiViii 
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CALCULATION  OF  THE  JACOBIAN  MATRIX 

To  complete  the  differentiation  of  the  Hamiltonian  with  respect  to  the  mean  variables, 
we  need  the  Jacobian  matrix 

J(g,  ff,  7,  S',  a) 

J(tj,  £,  r,  a,  L) 

for  the  three  cases:  zonal  harmonics,  lunar  attraction,  and  solar  attraction.  It  is  clear  from 
THE  HAMILTONIAN  IN  TERMS  OF  MEAN  VARIABLES  section  that  we  can  subsume  the 
three  into  one  by  writing 


a  =  or  x'  +  or  y'  +  a  z' 

x  y J  z 

(3  =  Bx'  +  By'  +  (3  z' 

x  y  *■ 

7  =  +  7yy'  +  yzz' 


where  the  unit  vector 


( 


Am\ 

for  the  zonals,/  y'  I 

I  M  I 

\M  / 


A 

for  the  moon,  and |  y's  J  for  the  sun 

VI 


(and  constant  in  all  three  cases),  and  J,  a,  ax,  ay,  a.z,  0x,  j3y,  (3z,  yx,  7y,  and  yz  are  the 
same  for  all  three  cases. 

If  we  define 

vi  =  -  T?o)/2Lf 


v2  =  -  (^o  +  tjt)/2LJ 


and  remember  that  c,  =  V(l+f)/2L,  c2  =  V (l+p)/2Lf 
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APPENDIX  G 

CALCULATION  OF  THE  JACOBIAN  MATRIX 
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Thus. 

C'm  (a.  0)  =  a  Cm_i  (a.  0)  -  0Sm_f  (a.  0) 

Sm  (a-  P)  =  a  Sm_)  (a.  j3)  +  /3Cm_(  (a.  0) 
and  the  starling  values  are 

C()  la.  0)  =  1  S0  (a.  0)  =  0 

Differentiation  from  the  definition  is 

ac 

r~~  («•  j3)  =  m  C  (a.  0) 
da  1 

ac 

(a-  <3)  =  -  m  Sm_,  (a.  (3) 
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Derivation  of  the  results  below  is  parallel  to  that  in  the  section  above. 


A™+)  (x)  =  (-)" 


(n-m)  !  Pn  <x> 


n!  x  e 


Upon  substitution  for  Pm  (see  the  previous  section), 


(n-m)  ! 

A™+1  (x)  =  (-)m  -  — —  x“n  +  m  Qnm  (x) 


n! 


The  recursion  relation  is 


n(n+l)  A™+2  (x)  =  n(2n  +  1)  A™+1  (x)  -  -+— 2  — A“  (x) 


with  starting  values 


A™  (x)  =  0 


A™  +  l  <*>  =  (") 


m  Qmm  M 


m! 


and  starting  recursion 


2m-  1 


A™>1  <X>  =  ~  —  Am  ” 1  A?  W  =  1 


Differentiation,  from  the  function  definition,  is 

A™+1  =  "  ~~x~  (An'+ 1  (X)  +  An'++ 1‘  <X» 


RECURSION  ON  C  (a,  0)  AND  S  (a,  0) 


These  coefficients  are  defined  by 


(a  +  i0)m  =  Cm  (a,  0)  +  i  Sm  (a,  0) 


so  that 


Cm  (a,  0)  +  i  Sm  (a,  0)  =  (a  +  i  /3)(Cm _ ,  (a,  0)  +  Sm_,  (a,  0)) 
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RECURSION  ON  Bm  (x) 

n 

These  are  almost  the  coefficients  defined  by  Cefola  and  Broucke  (Reference  F 
68-80).  but  multiplied  by  x2n+1,  which  simplifies  the  calculations  required. 


B 


m 

n  +  1 


(X) 


n! 

(n+m)  ! 


(x) 


(eq.  76) 


(Note  that  this  Pm  has  a  difference:  it  is  as  defined  in  Reference  F-2: 

n 


P™  (x)  =  (x2  -  l)m'2  Q™  (X),  not  (1  -  x2)m/2. 


(x) 


n! 

—  .  +  m  +  1 

(n+m)! 


Q™  (x) 


Now  the  recursion  follows  from  that  for  Qm  above. 

n 


(n-m+1)  (n+m+1)  B™+2  (x)  =  (n+1)  x2  l(2n+l)  B™+1  (x)  -  nB™  (x)] 
with  starting  values  B™  (x)  =  0.  B™  +  (  (x)  =  x2m  +  1/2m. 
and  the  starting  recursion  B™  +  )(x)  =  (x2/2)  B™ " 1  (x)  ,  B°  (x)  =  x 


Differentiation  of  the  function,  from  its  definition,  is 
d  n+m+ 1 

-  B™  ,  (x,  =  ——  (B-tl  (x)  +  Bnm;,>  (x)) 

RECURSION  ON  Am  (x) 

n 

These  are  the  coefficients  described  by  Cefola  and  Broucke  (Reference  F- 
54-67).  but  with  an  additional  factor  of  (-  l)m. 


1  Haul  J.  Cefola  ami  Roger  Broucke,  AIAA  Paper  No,  75-9  (no  date) 

2Milton  Ahramowitz  and  Irene  A,  Stegun,  Eds.,  “Handbook  of  Mathematical  H  unctions''  (Oover). 
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RECURSION  ON  Q™  (7) 


The  function  is  defined  by  Q™  (7)  = 


P  (7)  with  P  the  Legendre  polynomial  of  the 


n  d7m  n  n 

first  kind.  The  recursion  relationship  is  easily  obtained  by  differentiating  the  recursion  relation¬ 
ship  for  the  Legendre  polynomial: 


(n  +  2  -  m)  Qn+™  (7)  =  (2n  +  3)  7  Qn+"J  (7)  -  (n+m+I)Q™  (7) 
with  starting  values  Q™  +  1  (7)  =  0 

Q”  (7)  =  (2m  -  1)!! 

Qm+";  (7)  =  7  (2m  +  1)!! 


where  (2m  +  1)!!  =  1.3. 5. 7  .  .  .  .(2m  +  I).  The  starting  values  may  themselves  be  determined 
recursively: 

Q™  (7)  =  (2m  -  1)  Qm  _  t  m  " 1  (7) 

Qm+"J  (7)  =  (2m  +  1)7  Q™  (7) 

Q°  (7)  =  1 

Differentiation,  from  the  function  definition,  is 
±  Qn,  (7)  =  Q^1  +  (7) 


RECURSION  OF  D 

n 

(n-m)  ! 

Dm  is  defined  by  - -  Qm  (0)  and  its  recursive  relationships 

n  (n+m)  !  " 

,  n-  m+ 1 

follow  from  those  of  Qm  :  Dm  +  2  =  -  - Dm 

n  n  n+m+2  n 

Starting  values  are  Dm  =  — - ,  D  ™  ,  =0 

°  m  2m  m!  m  + 1 


and  the  former  allows  the  recursion 


Dm 

m 


—  Dm-[ ,  with  D°  =  1 
2m  m'‘  0 
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THE  SECOND-ORDER  SHORT-PERIOD  CORRECTION  TO  L 

The  last  equation  in  the  MEAN  ELEMENTS:  THE  VON  ZEIPEL  TRANSFORMATION 
shows  the  transformation  from  L  to  L'  with  the  second-order  correction.  The  portions  of  this 
involving  S2  are  calculated  in  the  subroutine  ELTRAN,  lines  273-278.  The  term  remaining, 
(L3/p2)(F2  -  F*),  is  more  conveniently  calculated  in  the  subroutine  DERIV,  lines  621-643, 
where  F*  is  already  available.  The  calculation  of  F2,  the  instantaneous  second-order  portion  of 
the  Hamiltonian,  is  as  follows: 

(1)  The  zonal  harmonics  contribution  to  F2. 

(AF),Z  =  -  E  (p/ r)  Jn  (R /r)n  P„  (sin  p) 

where  r  =  L2  f2  /pZ  and 

sin  (3  =  sin  i  sin  (f  +  g)  =  c2  (a  sin  u  +  t  cos  u) 

(2)  The  lunar  contribution  to  F2. 

(AF)1M  -  £  <«/rM)  <7r„)”  P.  (sin  ?') 

n  —  i 

where  r  is  as  above  and 
sin  p'  =  x'M  (cos  u  -  c2  r  sin  p/(  1  +  p)> 

+  yj^  (sin  u  -  c2  a  sin  j3/(l  +  p)  +  z'M  sin  P 


The  distance  to  the  moon  is  rM ,  and  the  components  of  the  unit  vector  to  the  moon  xj^ 

/v  A  A  , 

y'  ,  z'  ,  in  the  coordinate  frame  L,  L2  h,  and  c2  =  v(l  +  p/Lf. 

M  M 

(3)  The  solar  contribution  to  F2  is  exactly  the  same,  except  that  the  solar  coordinates 
are  substituted  for  lunar  coordinates. 

Since  F2  is  second-order,  it  makes  no  difference  whether  it  is  calculated  in  terms  of  mean 
or  osculating  variables. 

For  a  derivation  of  the  formulas  above  see  THE  HAMILTONIAN  IN  TERMS  OF  MEAN 
VARIABLES  section. 
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dS,  =  7a  |  U  j  +  p2^  dUfi  +  ((1  +  p)/4)  (o2  -  r2)dSup  + 
+  Id  +  2  p  -  5  p2)  Ue  + 

((3  +  5  p)/4L  ?)  ((a2  -  r1)  Sup  +  2  ar  Cup)]  (dL  -  £ 
-  (2  p  U8  +  (l/4Lf)  ((a2  -  r2)  Sup  +  2 ar  Cup))  (a  da  + 
+  (0  +  p)/2)  ((a  da  -  t  dr)  Sup  +  (a  dr  +  rda)  Cup)\ 


2  ar  dC  ) 

up' 


d£  -  77  dr?) 

r  dr) 


The  transformation  between  mean  and  osculating  coordinates  is  taken  from  these  equations 
and  is  contained  in  the  subroutine  ELTRAN,  lines  215-271. 
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In  exactly  similar  fashion. 


dS 

u  r 


(Z  cos 


2u/LJ3 )  l(Z  +  1) 


(1  +  n 


(sin  u  d£  +  cos  u  dr?) 


1  + 

1  +  ? 


e  sin  f(£  d£  +  r?  drj))  +  (1  +  f  +  s2)  (-tj  d£  +  £  dr?)  ] 


+  c  [sin  u  +  (1/3)  sin  3  u)d£  -  (cos  u  -  (1/3)  cos  3  u)dr?| 

Su  -  sin  2  u  (jL 

-  — -  (£  d£  +  77  dri)  +  - [-  2Z(Z  +  1)  cos  2  u.e  sin  f 

2L(1  +  £)  L?(l  +  r> 

-  (Sup  -  sin  2  u)[  +  2Z3  cos  2  u  dX/£3 


There  is  more  differentiation  to  be  done.  We  may  write 


S,  =  t'2  Lf  [(-  j  +  p2)  Ut.  +  ~~~  "  r2>  s„„  +  2  or  C„J| 


up 


u  p 


7:  (Lf)4 


i-Ly  +i  -L_ 

3  (Lf)3  <Lf)5  /  1  4  \(Lf)4 


<Ln5) 


((a2  -  r2 )  S  +  2  or  C  ) 


u  r 


U  p 


and  observe  that  7'2  (L  f)4  is  a  constant  and  that 


Lf  =  L.  -  (£2  +  i?2  )/2 

pL£  =  Lf  -  (a2  +  r2 )/ 2 


so  that  differentiating  S,  with  respect  to  L.f  and  p Lf.  and  Lf  and  pl_£  with  respect  to  th 
variables  yields  the  required  derivatives. 


END 
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